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1.  Introduction 


a 


The  research  covered  by  this  report  is  part  of  an  effort  to  develop 
techniques  for  the  evaluation  of  interpolation  methods  in  Digital  Elevation 
Models  (DEM*s). 

In  order  to  evaluate  the  fidelity  of  an  interpolation  method,  one  needs  to 
know  true  heights  in  an  area.  These  true  heights  can  then  serve  as  the  basis 
for  a  comparison  with  height  obtained  by  some  interpolation  method. 

In  practice,  one  deals  with  varying  types  of  terrain.  To  try  to  simulate 
conditions  encountered  in  practice,  one  therefore  needs  true  heights  for  a 
group  of  surfaces.  These  surfaces  should  represent  a  selection  of  the  types  of 
surfaces  which  may  be  encountered  in  practice. 

In  order  to  obtain  a  true  height  at  any  point  on  a  surface,  we  need  a 
synthetic  surface.  That  is,  a  surface  defined  by  a  mathematical  equation  where 
height  (or  elevation)  is  a  function  of  position  (x  and  y  coordinates). 

This  report  deals  with  methods  for  the  evaluation  of  contour-to-grid 
interpolation  methods  in  a  DEM.  It  is  in  essence  a  feasibility  analysis  and 
should  be  considered  as  an  effort  to  develop  an  initial  plan  for  the  evaluation 
of  DEM  interpolation  techniques. 

This  research  was  done  in  the  framework  of  the  Scientific  Services 
Program  -  STAS  and  the  specified  objectives  and  specific  tasks  are  given 
below. 

Objective: 

The  objective  of  this  short  term  analysis  is  to  develop  evaluation  tools  and 
methodqlogies  to  evaluate  contour-to-grid  interpolation  algorithms. 

Specific  Tasks: 

a.  Review  literature  regarding  relevant  research  in  the  area  of 
contour-to-grid  interpolation  evaluation  using  analytical  approximations. 

b.  Evaluate  techniques  for  generating  synthetic  surfaces  to  approximate 
Digital  Terrain  Elevation  Data  (DTED). 

c.  Develop  numerical  and  non-numerical  techniques  for  evaluating  the 
accuracy  of  contour-to-grid  interpolation  algorithms. 

d.  Demonstrate  the  methodology  for  creating  synthetic  surfaces  approximating 

DTED  and  evaluating  the  interpolation  techniques.  This  will  be  done  on  a 
small  sample  data  set.  — _ 
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Itelevant  literature  waa  reviewed  in  Journals»  conference  proceedingst 
reporta,  etc.  Referencea  related  to  the  reaearch  topic  were  aelected  and  are 
given  in  Appendix  A.  Each  reference  ia  accompanied  by  n  brief  deacription  of 
the  main  topiea  dealt  with. 

Nothing  apecific  waa  found  in  the  literature  dealing  with  interpolation 
methoda  in  a  OEM  where  the  OEM  aiay  repreaent  varying  typaa  of  terrain. 
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3.1  What  ia  a  Synthetic  Surface? 


Unlike  original  topographic  aurfacea,  aynthetic  aurfacea  are  thoae  aurfacee, 
that  can  be  fully  deacribed  by  an  analytical  function  of  the  form: 


z  =  f(x.  y) 


(3-1) 


Depending  on  the  mathematical  formulation,  aynthetic  aurfacea  vary  in 
complexity.  Such  eqxjationa  may  be  in  the  form  of  a  seriea  of  polynomiala: 


z  =  ax  +  by  +  cx*  dy*  exy 


(3-2) 


or  a  aeriea  of  trigonometric  functions: 


z  =  sin 


2n  ix 


(3-3) 


with  coainea  inatead  of  ainea  or  combinations  of  them.  On  the  other  hand 
covariance  functions  can  be  explicitly  formulated  (Loon,  1982)  and  from  these 
synthetic  surfaces  can  be  generated. 

All  these  methodologies  are  not  unknown  in  Cartography.  Morrison  (1971) 
has  generated  four  synthetic  surfaces  for  cartographic  use. 

The  way  in  which  the  aynthetic  surfaces  are  generated,  provides  them  with 
some  characteristics  which  do  not  alwaya  exist  in  topographic  surfaces.  They 
may  be  smooth,  continuous  (no  breaklines),  or  periodic  and  always  free  of 
measurement  errors.  This  makes  the  synthetic  surface  easier  to  handle  than 
an  original  topographic  surface. 


ismm. 


The  usefulneee  of  synthetic  surfacesy  however,  lies  on  the  fact  that  since 
they  are  strictly  defined,  they  can  provide  absolute  measures.  For  example, 
the  analytical  functions  used  in  this  study  produce  surfaces  which  are  smooth 
and  continuous  and,  being  free  of  any  measurement  error,  can  provide  an 
absolute  error  of  point  accuracy,  simply  by  comparing  the  height  value  at  the 
interpolated  point  to  its  true  value  aa  obtained  by  the  analytical  function. 
This  big  advantage  makes  the  use  of  synthetic  surfaces  a  necessary  tool  in 
this  study. 


Given  a  DBM,  where  the  heights  of  the  points  are  given  in  a  rectangular 
grid  pattern,  a  way  to  construct  a  synthetic  surface,  that  simulates  the 
original  terrain,  is  a  Fourier  aeries  expansion  of  the  given  DEM.  The  Fourier 
series  expansion  is  used  for  generating  synthetic  surfaces  in  this  project 
because  it  is  useful  for  collecting  information  relating  to  the  charactistics  of 
the  surface.  Future  studies  will  indicate  how  useful  this  method  will  be  for 
definitive  synthetic  surfaces,  but  other  methods  for  generating  synthetic 
surfaces  could  be  used.  An  in-depth  evaluation  of  all  the  possible  techniques 
for  generating  synthetic  surfaces  ia  a  project  on  its  own.  The  time  constraints 
in  this  project  have  not  allowed  us  to  concentrate  exclusively  on  this  aspect. 
The  synthetic  surfaces  generated  for  this  project  by  the  Fourier  method  have 
been  adequate  for  demonstrating  the  evaluation  of  interpolation  techniques. 
This  is  a  method  widely  used  (in  Geology,  and  Gravity  field  studies,  for 
example)  and  is  also  demonstrated  in  Morrison's  surface  III.  (Morrison  1971). 

The  basic  assumptions  underlying  the  Fourier  expansion  of  a  function  f(t), 
are  given  by  the  Dirichlet  conditions  (M.  Bath,  1974): 

1.  f(t)  should  be  periodic,  i.e.  f(t)  =  f(t  *  2n),  where  2n  is  the  period.  If  f(t) 
is  not  periodic,  but  defined  over  a  finite  range,  the  sum  of  the  sinusoidal 
terms  will  still  converge  to  f(t)  over  the  defined  range.  Outside  the  range, 
the  sum  will  represent  repetition  of  f(t). 

2.  f(t)  should  be  at  least  sectionally  continuous,  vdth  at  most  a  finite  number 
of  discontinuities  and  finite  jumps. 

3.  f(t)  should  possess  a  finite  number  of  maxima  and  minima. 

n 

4.  The  integral  J  f(t)  dt  should  be  convergent. 

-n 

According  to  Fourier's  theorem,  any  function  f(t)  satisfying  these 
conditions,  can  be  expressed  as  a  sum  of  a  finite  number  of  sinusoidal  terms. 
However,  it  must  be  noted  that  the  DiHchlet  conditions  are  sufficient,  but  not 
all  of  them  are  necessary.  This  means  that  the  Fourier  theorem  is  valid  for  a 
much  wider  class  of  functions. 


In  the  two-dimensional  case  of  the  Fourier  expansion  -  in  which  we  are 
interested  in  Cartography  -  a  function  z(x,  y)  can  be  expressed  as  follows 
(Davis,  1973): 


Zij(x.  y)  = 
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(3-4) 

where  we  assume  discrete  (n  x  m)  data  and  X,,  X,  refer  to  the  fundamental 
wavelengths  along  the  two  mutual  perpendicular  axes  x  and  y. 

Another  useful  representation  of  the  Fourier  series  is  in  complex  form. 
Assuming  again  two-dimensional  discrete  data,  the  representation  of  the  Fourier 
series  in  the  data  domain  is: 


N-l 

Z  =  f(x,  y)  =  ^ 
n=0 


M-1 

}  F(n,  u) 
m=0 


i2n 

e 


(3-5) 


where:  N  =  number  of  samples  along  x-direction 

M  =  ntjmber  of  samples  along  y-direction 


and  the  complex  Fourier  coefficients  F(n,  m)  are  orthogoned. 

Another  important  notion  is  the  Fourier  Transform  defined  as: 


_i2Tr  [M  +  EX] 
I  N  mJ 


(3-6) 


F(n,  ■) 


N-1 

n=0 


M-1 

YZ 

■=0 


The  term  spectrum 


is  often  used  as  a  synonym  for  the  Fourier  Transform. 


Using  Fourier  analysis  methods,  we  are  transforming  data  from  the  space 
domain  (X,  Y,  Z  coordinates)  to  the  frequency  domain  (ranking  of  frequencies 
present  in  the  surface).  In  the  space  domain,  we  can  calculate  tk.e  auto¬ 
correlation  (or  covariance)  function  of  points  in  a  DEM.  This  function  tells  us 
something  about  the  correlation  of  the  heights. 


UJ 

o 

u 


UJ 

a: 

a: 

o 

LJ 


L  n  G 


o 

u 


c 

1C 

o 


LOG 


Figure  1 


Autocorrelation  Curves 


Figure  2 


For  example,  in  Figure  1,  the  autocorrelation  curve  indicates  that  the 
heights  of  adjacent  points  tend  to  be  the  same  -  that  is,  we  have  a  fairly 
smooth  terrain.  On  the  other  hand,  in  Figure  2,  the  autocorrelation  curve 
indicates  that  height  of  adjacent  points  differ  greatly  -  that  is,  they  have  a 
low  correlation  and  therefore  the  surface  described  by  this  curve  is  rough  or 
complex.  The  autocorrelation  curve  deals  with  the  space  domain  and  the 
equivalent  concept  in  the  frequency  domain  is  the  power  spectrum  or  the 
power  density.  Power  spectra  are  frequently  used  in  geophysics  and  they  may 
be  calculated  in  two  ways: 


1.  First  calculate  the  autocorrelation  function  and  then  take  the  Fourier 
transform  of  this  function. 


2.  Calculate  the  Fourier  transform  of  the  original  surface  and  then  take  the 
square  of  its  absolute  value. 

The  power  spectrum  gives  us  information  about  the  contribution  of  the 
high  and  low  frequency  components  as  a  function  of,  in  our  case,  the  elevation. 
Surfaces  with  strong  low  frequencies  and  relatively  unimportant  high 
frequencies  tend  to  be  smooth  with  small  local  fluctuation.  On  the  other  hand, 

5 


surfaces  whose  high  frequencies  contribute  most  of  the  total  power  tend  to 
have  extreme  local  fluctuations. 


Other  researchers  (Frederiksent  et  al.,  1984)  have  used  the  (2-D)  power 
spectrum  of  profiles  for  terrain  classification.  There  may  be  much  scope  for 
using  the  (3-D)  pqwer  spectrum  of  surfaces  for  the  classification  or  matching 
of  terrain  regions.  This  is  an  area  that  still  needs  investigating. 

Therefore  from  the  Fourier  transformi  the  power  densities  are:  (Engelis, 
1983) 


IF(0,  0)1* 

:  for  n=B=0 

4IF(n,  m)l* 

:  for  n=l,  . . .  Li-1. 

n— 1 ,  ...  La— 1 

4IF(L»,  LjI* 

:  for  n=Li,  ■=L2, 

N,  M=odd 

2IF(L,,  Lai* 

for  n=Li ,  B=La 

N,  M=even 

(3-7) 

The  total  power  is  given  by: 


P 


L. 


r 


(3-8) 


Two  considerations  at  this  point  are  very  important: 

1.  Depending  on  the  data  sampling,  the  highest  recoverable  frequency  in  each 
direction,  is  called  the  Nyquist  frequency  (denoted  by  L,  for  x  and  L2  for 
y  axis)  given  by: 


N±1  . 
2  • 

for 

N:odd 

N 

2  • 

for 

N:even 

(3-9) 


2.  "At  least  theoretically,  the  mean  value  of  the  original  data  is  not  siffected 
by  the  filtering  operation.  Indeed,  since  the  filter  coefficients  add  up  to 
unity,  and  in  theory  the  sampled  data  are  of  infinite  extent,  the  mean 
value  of  the  input  data  is  the  same  as  the  one  of  the  output  results.  This 
is  not  the  case  though,  for  the  actual  filtering  of  finite  samples,  where  the 
mean  is  also  subject  to  the  occuring  distortions.  So,  in  order  to  reduce  as 
much  as  possible  the  effect  of  these  distortions,  the  mean  value  of  the  data 
should  be  subtracted  before  the  filtering  process  and  then  added  back  to 
the  filtered  values."  (Engelis,  1983). 


Sinc«  th«  Pouriar  tranaform  ia  computed  for  the  DEM  using  the  Fast 
Fourier  Transform  Method,  the  power  spectrum  and  the  percentile  contribution 
of  each  frequency  component  to  the  overall  power  can  be  obtained.  A  kind  of 
low-paaa  filtering  can  then  be  applied,  only  the  significant  frequencies  are 
retained.  Their  significance  is  actually  tested  against  a  tolerance  value,  as 
describ'd  below. 


Assuming  that  Z(z,  y)  is  a  sequence  of  normal,  independent  (m, 
variables,  the  Fourier  coefficients  a,  b,  c,  d  being  linear  combinations  of  Z(x, 
y)  will  be  norsutUy  distributed.  These  coefficients  are  also  independent,  since 
the  sine  and  cosine  functions  are  orthogonal  in  the  basic  interval. 


Assusm  that  iP(n,  m)l  ia  the  amplitude  of  the  frequency  f^,  m  with  power 
P(n,  m)  o**  variance  e^in.  The  significance  of  this  frequency  can  be  tested  with 

a  t-test  as: 


P{t  >  toF.a)  =  1  -  a 


(3-10) 


where 


^DF,a*  critical  value  as  obtained  from  t-distribution  tables 


a:  type  I  error 


OF:  degrees  of  freedom 


and  the  null  hypothesis  to  be  tested  is: 


Hq  =  P(n,n)  =  0 


with  the  alternative: 


F(n,.)  t  0 


If  (3-10)  holds  we  reject  Ho,  which  means  that  the  frequency  fn,m  should  be 
retained  as  significant  at  (l-«i)X  significance  level. 


An  important  point,  here,  would  be  the  determination  of  the  degrees  of 
freedom.  Each  frequency  corresponds  to  4  degrees  of  freedom  (vi=4,  in  the 
two-dimensional  case),  as  by  virtue  of  the  Fourier  theorem  there  are  four 
constants  (a,  b,  c,  d)  needed  to  characterize  each  frequency  component  in  a 
waveform.  (Bath,  1974).  On  the  other  hand,  when  spectra  are  calculated 
directly  by  the  Fast  Fourier  Transform  method,  each  raw  estimate  has  2 
degrees  of  freedom  (in  the  two-dimensional  case).  Averaging  over  N  estimates 
(when  taking  the  transforms  with  respect  to  one  of  the  variables)  yields  new 


estimates  with  2N  degrees  of  freedom.  Then,  the  transforms  are  taken  with 
respect  to  the  second  variable,  and  finally  we  have  2N  x  2M  (va=4NM)  degrees 
of  freedom. 


After  retaining  only  the  significant  frequences,  the  Inverse  Fourier 
Transform  can  be  used  to  obtain  the  filtered  surface  in  the  data  domain.  Thm 
surface  is  now  strictly  known,  since  the  complex  Fourier  coefficients  are 
determined,  and  can  serve  as  a  synthetic  surface 


3.3  Generation  of  Synthetic  Surfaces  Used  in  This  Study 

The  synthetic  surfaces  used  in  this  study  have  been  generated  using 
Fourier  series  expansion  of  six  original  topographic  surfaces.  As  explained 
previously,  the  density  of  the  data  prohibit  the  recovery  of  terms  with 
frequency  higher  than  a  boundary  value  (Nyquist  frequency).  This  means  that 
the  Fourier  aeries  are  not  of  infinite  length.  Moreover,  not  all  the  frequencies 
are  significant  and  contribute  the  same  way  in  the  given  surface.  Therefore 
performing  a  t-test  of  significance,  we  can  furthermore  truncate  the  Fourier 
series  retaining  only  the  significant  terms.  This  can  provide  us  with  the 
formulation  of  a  Fourier  series  in  exponential  form  which  keeps  the 
characteristic  of  the  original  terrain  and,  mathematicadly  generates  a  synthetic 
surface. 

For  this  project,  we  chose  to  generate  synthetic  surfaces  which  can  be 
described  by  a  Fourier  series  up  to  the  third  harmonic  surface.  The  third 
harmonic  surface  gives  49  terms  in  the  equation  for  the  synthetic  surface. 
Higher  order  harmonic  surfaces  would  greatly  increase  the  complex 
mathematical  manipulation.  Time  did  not  allow  us  to  include  these  higher  order 
harmonic  surfaces.  Spectrum  analysis  has  shown  that,  in  the  majority  of  cases, 
most  of  the  power  of  the  original  surface  is  contained  in  the  harmonic  terms  up 
to  the  third  harmonic  surface.  Comparing  the  power  in  the  first  three 
harmonic  surfaces  with  the  power  in  the  first  nine  harmonic  surfaces  in  Table 
1  will  confirm  the  previous  statement. 

The  synthetic  surfaces,  as  can  be  seen  in  the  figures  below,  do  not  look 
like  "terrain  types"  which  one  would  expect  in  mapped  topography.  If  the  aim 
of  this  project  was  to  classify  terrain  types,  then  this  would  be  a  serious 
fault.  In  terrain  type  classification  we  would  expect  to  find  areas  with 
breaklines  and  varying  degrees  of  surface  complexity.  In  this  project, 
however,  we  restrict  ourselves  to  "interpolation  regions”,  that  is,  regions 
where  interpolation  algorithms  will  be  applied.  The  boundaries  of  these 
interpolation  regions  are  the  breaklines  because  a  basic  assumption  in 
interpolation  (especially  in  cartography)  is  that,  one  never  crosses  a  breakline. 
Therefore,  the  fact  that  the  interpolation  region  does  not  look  like  terrain  or 
the  fact  that  the  interpolation  region  show  some  periodic  pattern  will  not 
invalidate  our  suggested  evaluation  methods.  In  Clarke  (1982)  it  is  shown  that 
when  one  applies  an  interpolation  algorithm  to  topographic  surfaces  the  results 
are  better  than  what  one  gets  from  synthetic  surfaces.  This  could  possibly  be 
due  to  the  fact  that  topographic  surfaces  carry  more  information  in  the  form  of 
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additional  height  values  at  critical  points. 

Throughout  this  project  one  needs  to  bear  in  mind  that  interpolation 
regions  are  being  discussed  and  not  terrain  types. 

The  first  five  surfaces  have  been  based  on  the  15*  topographic  series, 
scale  1:62,500,  and  cover  areas  of  9.5  km  z  9.5  km.  A  DEM  has  been 
constructed  (by  estimation)  from  countours  for  each  one  of  them,  consisting  of 
20  X  20  points  500  m  apart. 

The  sixth  surface  has  been  based  on  the  7X*  topographic  series,  scale 
1:24,000  and  covers  an  area  of  4.75  x  4.75  km.  The  DEM  which  has  been 
constructed  consists  of  20  x  20  points  250  m  apart. 

The  surfaces  based  on  the  above  maps  are  given  below.  The  name  of  the 
sheet  is  given  together  with  the  U.T.M.  coordinates  of  the  South-West  corner  of 
the  square  from  which  the  elevations  were  taken. 

These  synthetic  surfaces  do  not  resemble  the  original  surfeu:es  (as 
represented  on  the  topographic  maps).  The  reasons  are  given  above  and  we 
need  to  add  that  relatively  few  points  were  used  to  sample  the  original 
surfaces  -  points  500  m  apart  for  the  first,  five  surfaces  and  250  m  apart  for 
the  sixth  surface.  This  means  that  wavelengths  smaller  than  1000  m  (twice  the 
sampling  distance)  in  the  first  five  surfaces  and  smaller  than  500  m  in  the 
sixth  surface  cannot  be  represented.  The  topographic  map  surfaces  were  used 
as  a  basis  for  generating  synthetic  surfaces  of  varying  complexities.  This 
complexity  is  defined  by  the  power  spectra  of  the  surfaces  (see  Table  1). 

SURFACE  1:  Mineral  King,  CALIF.,  Qudrangle 
x  =  355,000  ■  E 
.  y  -  4,026,560  a  N 

SURFACE  2:  Georgetown,  COLO.,  Quadrangle 
X  =  447,550  a  E 
y  =  4,374,800  a  N 

SURFACE  3:  Sundance,  WYO.,  Quadrangle 
X  =  540,000  a  E 
y  =  4,, 899, 600  a  N 

SURFACE  4:  Reedsport,  OREG. ,  Quadrangle 
X  =  409,800  a  E 
y  =  4,816,900  a  N 

SURFACE  5:  New  Haven,  ILL-IND-KY,  Quadrangle 
X  =  399,000  a  E 
y  =  4,188,800  a  N 

SURFACE  6:  Horse  Cave,  KY. ,  Quadrangle 
X  =  589,000  a  E 


The  synthetic  surfaces  (that  is,  as  defined  by  an  equation)  are  numbered  1 
to  6.  The  "surfaces"  as  defined  by  the  20x20  array  of  sampled  points  (as 
described  above)  are  numbered  lA  to  6A. 

The  program  PPTANAL  (see  Appendix  B)  was  used  to  transform  the  data  to 
the  frequency  domain  and  perform  a  spectrum  analysis. 

Table  1  shows  the  power  spectrum  of  each  DEM  and  the  contribution  of 
each  harmonic  surface  to  the  total  power. 


The  six  surfaces  used  are  of  differing  compleziiies.  This  can  be  seen  bjr 
visual  insiMction  (that  ist  do  the  contourSt  in  a  general  sense,  show  complex  or 
simple  form)  or  by  analysing  their  spectra.  A  general  ranking  of  the  order  of 
complexity  is  as  follows: 

Surface  ZA 
Surface  lA 
Surface  4A 
Surface  3A 
Surface  6A 
Surface  5A 


The  complexity  will  be  defined  by  the  power  spectrum.  An  examination  of 
Table  1  shows  that  the  above  ranking  holds  true  when  one  examines  the  first 
three  harmonic  surfaces  or  when  one  examines  the  first  nine  harmonic 
surfaces. 

As  we  will  note  below,  the  complexity  of  the  surface  will  affect  the 
accuracy  of  the  interpolation. 


3.4  Generation  of  Synthetic  Surfaces  up  to  3***^  Harmonic  Surface 

As  described  in  previously,  a  Fourier  series  expansion  up  to  the  third 
harmonic  surface  follows  equation  (3-4)  where  n=m:0,  3,  and  this  gives  49 
non-zero  harmonic  coefficients.  These  unknown  coefficients  have  been  found 
for  each  surface  using  the  least  squares  adjustment  program  FOURIER.  (See 
Appendix  B).  The  equations  of  all  surfaces  were  obtained  and  are  included  in 
Appendix  C. 

The  contouring  algorithm  used  for  all  the  output  in  this  project  is  the 
Geodetic  Science  Plotting  Package  (GSPP)  as  described  in  Sunkel  (1980). 

The  following  figures  show  the  power  spectra  of  the  six  surfaces.  The  3-D 
representations  show  only  one  quarter  of  each  spectrum  -  the  other  three 
quarters  are  symmetric. 

In  some  of  the  related  references  (for  example  Bath  1974,  Jacobi  1980)  the 
spectrum  is  given  in  graphs  of  log  P  vs  log  1/X,  where  P  is  the  power  of  each 
frequency,  t,  and  X  is  the  corresponding  wavelength.  This  is  siainly  used  to 
find  the  sIoim  of  the  spectrum  which  is  a  characteristic  of  different  terrain 
types.  The  following  logarithmic  scales  show  the  plots  of  log  P  vs  f.  This 
exaggerates  the  high  frequencies  which  normally  have  low  amplitudes. 


An  euaination  of  the  following  fitfurea  ehowe  that  when  the  power 
epectrum  ie  plotted  at  a  logarithmic  ecale  then  more  details  of  the  spectrum 
becosie  visible.  For  future  research  this  could  make  it  easier  to  classify 
terrain  or  interpolation  regions.  It  may  also  be  easier  to  match  existing 
terrain  spectra  with  previously  defined  criteria  spectra.  For  this  purose  one 
could  fit  a  plane  to  the  spectrum  and  examine  the  slope  and  direction  of  this 
plane.  For  the  2-D  casei  logarithmic  scales  were  used  by  Frederiksen  et  al. 
(1984). 

The  following  figures  of  power  spectra  give  us  general  impressions  of  the 
character  of  the  interpolation  regions.  For  example,  comparing  figures  7  and 
11  we  can  note  that  the  terrain  represented  by  figure  11  is  smoother,  less 
rugged  and  more  homogeneous  than  the  terrain  represented  by  figure  7.  The 
homogeneity  of  surface  3  is  seen  from  figure  11  where  the  power  spectrum 
shows  that  in  both  the  X  and  Y  directions  the  same  frequency  components  have 
approximately  the  same  power. 

It  is  clear  to  us  that  much  work  needs  to  be  done  in  the  area  of 
classification  of  surfaces.  For  this  project,  we  will  classify  the  complexity  of 
the  synthetic  surfaces  according  to  the  value  obtained  for  the  power  of  the 
harmonic  surface.  Ordering  of  complexity  according  to  the  elements  of  the 
spectrum  (as  shown  in  the  following  figures)  still  needs  investigation.  Fitting 
a  plane  to  the  3-D  spectrum  may  be  one  way  of  accomplishing  this  ranking  or 
ordering. 


Noraal  Scale  Logarithaic  Scale 


In  order  to  generate  a  raster  image  of  synthetic  surface  contours,  a  blank 
raster  matrix  of  size  101x101  points  25  m  apart  is  set  up.  Then  for  each  pixel 
(x,  y)i  the  equation  of  the  synthetic  surface  is  solved,  obtaining  a  zi  value  for 
the  height.  Then  setting  a  contour  interval,  every  z  value  close  to  a  contour 
value,  within  a  tolerance  limit,  takes  the  value  of  the  contour,  the  other  pixels 
are  zeroed  out,  and  thus  a  raster  image  is  produced.  For  example,  see  figure 
19(a). 

Some  errors  produced  by  this  procedure  will  be  added  to  the  errors 
produced  by  the  interpolation  algorithms.  Such  errors  are: 

1.  In  very  steep  areas  gaps  may  occur  in  the  contour  lines. 

2.  The  tolerance  values  used  have  already  introduced  an  error  in  contour  line 
determination. 

However,  no  attempt  has  been  made  to  correct  these  errors  in  contour  line 
determination.  The  reasons  are: 

1.  They  reflect  actual  situations  where  a  scanner  is  used  to  scan  an  existing 
contour  map.  Gaps  may  be  left  in  contour  lines  or  very  thick  contours 
may  be  produced,  depending  on  the  resolution. 

2.  Depending  on  the  algorithm  used  to  draw,  smooth,  or  generalize  the 
contours  on  a  contour  map,  errors  of  different  magnitudes  are  introduced. 


3.6  Complexity  of  the  yermonic  Surfaces 

As  mentioned  previously,  the  six  terrains  used  in  this  study  are  of 
different  complexity.  The  resulting  3^^  harmonic  surfaces  are  also  of  different 
complexity  as  shown  in  Table  3,  where  the  order  of  complexity  is:  SURFACE  2, 
SURFACE  1,  SURFACE  4,  SURFACE  3,  SURFACE  6  and  SURFACE  5  from  the  most 
complex  to  the  least  complex. 

Table  2  below  gives  the  mean  height,  variance  and  total  power  of  each 
surface  and  Table  3  below  shows  that  most  of  the  power  of  the  original  sampled 
terrain  is  included  in  these  synthetic  surfaces  when  they  are  expanded  to 
include  all  the  frequencies  existing  in  the  harmonic  surfaces  up  to  the  third 
order. 


Table  2 


I 


Surface 


Some  Surface  Attributes 
(the  array  of  sampled  points) 


Mean  Height 

Variance 

Total  Power 

10,372.023  ft. 

818,611.000 

1,127,147.00 

10,594.547 

995,733.437 

1,870,697.00 

4,885.422 

20,067.152 

19,860.60 

528. 100 

87,116.375 

132,066.87 

389.042 

360.502 

629.57 

765.375 

4,935.863 

5,795.30 

Table  3 

Power  of  the  Synthetic  Surfaces 

Contribution  of 

3*^^  h.s.  to  the  total  power 

power  of 

3*“'*  h.s. 

84» 

954,735.37 

89!t 

1,663,460.00 

80X 

15,842.46 

BIX 

81,075.62 

TJX 

483.34 

47» 

2,741.20 

The  complexity  of  these  surfaces  has  a  direct  effect  on  the  definition  of 
the  contour  lines,  as  will  be  discussed  below.  These  errors  are  therefore 
added  to  the  actual  errors  in  the  interi>olation  algorithms. 


4  Interpolation  Methods 

4.1  Interpolation  from  Contoura  to  Grid 

It  is  often  desired  to  estimate  the  value  of  a  function  in  locations  where  it 
has  not  been  sampled.  This  process  is  known  as  interpolation.  The  kind  of 
interpolation  to  be  used  must  be  such  that  takes  advantage  of  all  the 
information  contained  in  the  data,  while  at  the  same  time  does  not  produce 
artificial  results.  This  means  that  the  design  of  the  method  should  be  closely 
related  to  the  type  and  the  organization  of  the  input  data. 

Conversion  of  contour  lines,  in  digital  form,  to  a  rectangular  DEM  is  a 
problem  that  cannot  be  handled  by  some  general  interpolation  method.  The 
first  constraint  is  the  enormous  amount  of  data  obtained  during  digitization  or 
scanning  of  the  contour  lines.  A  study  of  the  published  contour-specific 
algorithms  suggests  that  the  most  accurate  results  should  be  obtained  if  the 
actual  interpolation  takes  place  along  the  line  of  the  steepest  slope.  In  the 
next  section  three  algorithms  are  presented  as  examples  and  an  effort  to 
evaluate  their  performance  and  efficiency  is  made. 

We  need  to  distinguish  between  "method"  and  "algorithm”.  A  method  is  a 
manner  of  procedure,  that  is,  a  description  of  how  something  could  be  done. 
An  algorithm  translates  a  method  into  a  set  of  well-defined  processes  which 
leads  to  the  solution  of  a  problem  in  a  finite  number  of  steps.  Algorithms  can 
be  turned  into  computer  programs  but  not  all  methods  can  be  turned  into 
algorithms. 


4.2  Contour  Specific  Algorithms 

The  problem  of  interpolation  seems  to  be  more  peculiar  in  the  case  of 
digitized  contours,  than  in  any  other  general  case.  There  are  a  number  of 
reasons  for  that: 

1.  This  field  has  not  been  explored  in  great  detail.  Algorithms  are  available 
in  the  literature,  (Roberts  1980,  Leberl  et  al.  1980,  for  example),  but  no 
evaluation  tools  are  reported  in  the  literature. 

2.  The  data  structure  of  digitized  contours  has  a  particular  form.  That  is, 
the  data  is  in  string  format  and  concentrated  along  lines  leaving  large 
areas  devoid  of  data. 

3.  Besides  the  height  of  the  points,  the  contours  give  additional  information, 
which  should  be  accomodated  in  the  algorithm.  The  contours  describe  the 
structure  of  the  surface. 

4.  Besides  the  contours,  other  lines  of  special  interest  (breaklines,  ridge 
lines)  and  spot  heights  provide  supplementary  information,  which  should 
not  be  ignored  in  the  interpolation  algorithm. 


In  this  study,  three  contour-specific  interpolation  algorithms  have  been 
evaluated  using  synthetic  surfaces.  In  the  next  sections  the  form  and  the 


structure  of  the  input  data  and  the  logic  of  interpolation  method  are 
explained. 


4.2.1  Input  Data 

Normally t  the  input  data  for  these  algorithms  consists  of  digitized  contours , 
breaklines  and  spot  heights.  After  the  digitization  in  vector  mode  all  the 
above  information  is  converted  into  a  raster  image  of  the  terrain.  In  this 
study,  since  synthetic  surfaces  have  been  used,  no  breaklines  or  spot  heights 
have  been  incorporated.  The  actual  input  data  consists  of  a  raster  image  of 
101x101  pixels,  describing  each  surface  by  contours  in  raster  mode.  The 
output  is  a  DEM  of  21x21  height  points,  125  m  apart  each  other  so  to  comply 
with  DTED  specifications  used  by  DMA.  (DMA  1979).  Other  parameters  should 
also  be  specified  in  these  algorithms,  i.e.  dimensions  of  raster  image  and  output 
DBM,  resolution  of  the  raster  (pixel  size  in  mm),  contour  interval,  and  ground 
resolution  of  the  pixels.  Finally,  at  least  approximate  height  values  for  the 
four  corners  of  the  output  DEM  should  be  provided  by  the  user.  For  more 
information,  refer  to  A.L.  Clarke,  (1982).  Figure  19(b)  illustrates  the  principles 
of  the  following  algorithms. 

4.2.2  LIXY:  Linear  Interpolation  Alon^r  X  and  Y  Axes 

This  method  implies  linear  interpolation  along  two  pre-specified  axes.  It  is 
computationally  efficient,  but  it  is  not  expected  to  be  very  accurate  especially 
in  the  case  of  steep  terrain.  The  algorithm  is  rather  simple:  the  height  at  the 
interpolation  point  is  found  by  linear  interpolation  between  the  points  on  the 
nearest  contours  in  the  X  and  Y  directions.  See  Appendix  B. 

The  whole  procedure  includes  the  following  steps: 

1.  Approximate  height  values  for  the  four  corners  of  the  derived  DEM  are 
provided  by  the  user. 

2.  The  interpolated  height  values,  for  the  points  of  the  DEM  along  the  four 
edges  of  the  map,  are  found  by  linear  interpolation  between  the  nearest 
contours  along  these  edges. 

3.  The  corner  and  edge  values  of  the  DEM  are  added  back  to  the  original 
contour  map,  as  already  gained  information. 

4.  The  corresponding  row  (x-direction)  and  column  (y-direction)  is  found  for 
the  next  interpolation  point. 

5.  If  the  point  happens  to  be  on  a  contour  line  this  contour  value  is  assigned 
to  the  point  as  its  elevation.  Otherwise  the  interpolation  value  is  found  as 
the  mean  of  the  two  values  resulting  by  linear  interpolation  between  the 
nearest  contours  along  the  x  and  y  axes. 


6.  Finally,  the  output  is  s  line  printer  map  of  the  raster  data  and  the  DEM 
point  locations,  and  a  listing  of  the  interpolated  DEM. 


4.2.3  LISS:  Linear  Iniervolation  Along  the  Direction  of  the  Steepest  Slope 

This  algorithm  is  a  variation  of  the  LIXY  algorithm  and  is  based  in  the 
first  place,  on  the  definition  of  the  direction  of  the  steepest  slope,  and 
secondly,  on  the  linear  interxiolation  between  the  two  nearest  contour  {mints 
along  this  profile. 

For  the  definition  of  the  direction  of  the  steepest  sloim  four  directions  are 
searched;  N-S,  W-B,  NW-SE,  NB-SW.  Then  the  direction  of  the  steepest  slope  is 
approximated  by  one  of  these  directions  that  show  the  steepest  sloim. 

The  steps  of  the  algorithm  are: 

1.  The  user  must  define  the  parameters  needed  in  the  program  and  provide 
approximate  height  values  for  the  four  corner  {mints  of  the  derived  DEM. 

2.  The  interpolated  values  of  the  DEM  {mints  along  the  edges  of  the  map  are 
found  by  linear  interpolation  between  the  nearest  contour  points  along 
these  edges. 

3.  The  corner  and  edge  points  of  the  derived  DEM  are  added  back  to  the 
raster  image. 

4.  For  every  DEM  point  in  the  interior,  eight  data  {mints  in  the  X,  Y  and 
diagonal  directions  are  located  and  the  profile  with  the  steepest  slo{>e  is 
found.  If  two  or  more  profiles  have  the  steepest  slope,  then  the  profile 
with  the  shortest  distance  is  chosen. 

5.  Along  this  profile  a  linear  interpolation  is  performed  and  the  interpolated 
height  value  for  the  DEM  is  computed. 

6.  The  program  output  consists  of  a  line  printer  map  of  the  raster  data  and 
the  DEM  point  locations,  a  listing  of  the  inter{mlation  directions  and 
heights  for  each  DEM  point,  and  a  listing  of  the  derived  DEM. 

This  method  differs  from  the  previous  one  in  that  the  inter{mlation  is 
performed  in  one  profile  only.  The  search  for  the  steepest  slope  makes  the 
LISS  algorithm  less  efficient,  in  terms  of  cost,  than  the  LIXY,  but  the  results 
obtained  by  this  method  are  expected  to  be  more  accurate. 

4.2.4  CISS:  Cubic  Interpolation  Along  the  Direction  of  the  Steepest  Slave 

This  is  the  most  complex  algorithm  since,  besides  the  search  for  the 
direction  of  the  steepest  slope,  it  requires  the  location  of  additional  {mints  and 
performs  interpolation  based  on  a  cubic  {mlynomial  or  spline  fitting.  The 
computational  time  is  also  greater  but  the  results  are  ex{}ected  to  be  even  more 
accurate  than  the  above  methods  especially  in  steep  terrains. 

The  proposed  model  here,  is  a  "moving  bicubic  spline  patch”  inter{mlation 
algorithm.  A  bicubic  spline  surface  is  local,  smooth  and  exactly  fits  the  input 
data.  S{mt  heights  (local  minima  and  maxima)  may  be  incorporated  by  setting 
the  first  derivative  of  the  surface  to  zero,  where  the  first  and  second 
derivatives  (slope  and  curvature)  can  be  manipulated  to  accomodate  breakline 


features.  Further  approximations  make  the  algorithm  more  efficient:  reduction 
of  the  bicubic  spline  surface  to  a  three-dimensional  cubic  space  curve  in  the 
direction  of  the  steepest  slope  and  finally  the  approximation  of  the  space  curve 
by  a  two  dimensional  profile  in  the  direction  of  the  steepest  slope  containing 
the  interpolation  point. 

The  "minimum  norm"  and  the  "best  approximation"  (Moritz,  1978)  properties 
of  the  cubic  spline  are  desired  in  the  DEM  interpolation.  The  end  conditions 
for  the  splines  are  essentially  arbitrary.  Natural  splines  assume  second 
derivatives  set  to  zero  at  the  end  points  which  is  equivalent  to  linearity  at  the 
extremities,  and  so  produces  flatter  curves.  Because  of  this  stability,  a  natural 
spline  with  conditions  is  used  by  the  CISS  algorithm.  Tests  have  also  shown 
(A.  Clarke,  1982)  that  "either  a  cubic  polynomial  or  natural  cubic  splines  may 
be  used  for  interpolation  in  the  central  segment,  but  only  the  spline  function 
is  suitable  in  an  end  segment.  The  similar  results  for  the  central  segment 
allow  the  choice  of  an  interpolation  function  for  that  segment  to  be  based  on 
computational  considerations". 

The  actual  steps  of  the  algorithm  are: 

1.  Parameters  needed  by  the  program  should  be  defined  and  approximate 
values  for  the  DEM  corners  must  be  provided  by  the  user. 

2.  For  the  edge  points  the  interpolated  height  value  is  found  using  cubic 
spline  interpolation  using  four  data  points  located  on  the  intersections  of 
the  map  edge  and  the  four  nearest  contours. 

3.  The  corner  and  edge  points  are  added  back  to  the  raster  data  to  provide  a 
boundary. 

4.  As  in  LISS  the  direction  of  the  steepest  slope  is  found. 

5.  Along  this  profile  the  interpolated  height  is  computed.  The  kind  of 
interpolation  used  depends  on  the  local  situation: 

a)  If  no  breaklines  of  edges  are  encountered  in  the  neighborhood,  a 
Newton  form  cubic  polynomial  using  the  four  nearest  contour  points  is 
fitted  to  the  data,  two  points  on  each  side. 

b)  If  a  breakline  of  data  edge  is  located  in  the  steepest  slope  profile, 
natural  cubic  spline  interi>olation  is  used. 

c)  If  both  ends  of  this  profile  are  on  a  breakline  or  data  edge,  then 
linear  interpolation  is  used. 

6.  The  interpolated  height  is  compared  to  the  adjacent  data  point  for  gross 
errors  detection. 

7.  The  program  output  consists  of  the  usual  line  printer  map,  a  listing  of  the 
derived  DEM  and  detailed  information  for  each  DEM  point  interpolation  (for 
example,  direction  of  steepest  slope,  raster  and  DEM  point  location,  terreun 
slope,  number  of  points  used  in  the  interpolation,  flags  for  blunders). 


5.  Accuracy  Techniques 

The  basic  scenario  is  as  follows.  Contours,  representing  a  surface,  are 
digitized.  An  interpolation  method  is  then  used  to  predict  points  on  a  regular 
grid.  As  we  are  using  synthetic  surfaces,  the  true  values  of  these  predicted 
points  are  known. 

By  comparing  the  predicted  values  and  the  true  values,  one  arrives  at  a 
number  of  measures  which  indicate  the  accuracy  of  the  prediction  over  the 
surface.  These  measures  are: 


mean  error; 

standard  deviation  of  error; 
maximum  absolute  error 

and  any  of  the  above  as  a  percentage  of  the  contour  interval. 

One  can  obtain  added  insight  into  the  accuracy  of  the  interpolation  method 
by  looking  at  a  map  of  the  residuals  (true  value  minus  predicted  vedue)  plotted 
in  their  correct  locations.  This  visual  aid  can  indicate  the  existence  of 
systematic  error,  as  well  as  other  isolated  errors  (or  perhaps  blunders)  in  any 
part  of  the  surface. 

To  understand  what  all  the  above  really  means  in  practical  terms,  we  have 
added  maps  showing  contours  produced  from  grid  values  which  were  predicted 
(by  the  interx)olation  method)  superimposed  on  contours  produced  from  true 
grid  values  (obtained  from  the  synthetic  surface  equation).  Ideally,  the  two 
sets  of  contours  should  be  plotted  in  different  colors  but  we  did  not  have  the 
facilities  to  do  this.  The  difference  between  the  two  sets  of  contours  gives  an 
impression  of  the  practical  effects  of  all  the  interpolation  errors. 


6.  Reaulta  Usin*  Small  Datasets 
6.1  Numerical  Evaluation  Results 


The  purpose  of  obtaining  a  DEM  from  digitized  contours  was  to  derive  a 
digital  form  of  a  surface  which  complies  with  the  accuracy  of  the  "Digital 
Terrain  Elevation  Data"  (DTED).  The  accuracy  of  the  DTED  eis  specified  by  the 
Defense  Mapping  Agency  (DMA  1979)  is:  "DTED  accuracy  is  specified  as  +12  m 
linear  error  relative  to  mean  sea  level  with  90X  confidence".  What  this  means 
statistically  is  that  the  probability  of  an  error  to  be  less  than  -12  m  and 
greater  than  +12  m  is  90%. 

Assuming  that  errors  follow  a  normal  distribution,  the  probability: 


P{xt  <  e  ^  Xj}  =  90X 


(6-1) 


yields: 

F(xi)  =  0.05 
F(xa)  =  0.95 


(6-2) 


where  F(xi)  is  the  cumulative  probability  density  function  of  the  normal 
distribution.  Using  the  tables  for  the  normal  distribution  we  get: 


X,  =  -1.6449.<r 


Xa  =  1.6449.<r 


(6-3) 


where  <r  =  the  standard  deviation  of  the  errors. 
Therefore  (6-1)  becomes: 


P{-1.6449  <r  *  e  *  1.6449  v}  =  90* 


(6-4) 


The  next  table  (Table  4)  gives  the  confidence  intervals  of  the  errors  for  the 
six  datasets,  using  the  three  interpolation  algorithms,  at  90X  significance  level. 


The  errors  mentioned  in  Tables  4  and  5  are  true  errors,  because  the 
synthetic  surface  equation  gives  the  true  values.  That  is,  the  error  is  the 
difference  between  the  true  value  (as  obtained  from  the  synthetic  surface 
equation)  and  the  interpolated  value  (as  obtained  from  the  particular  algorithm 
being  examined). 

We  need  to  bear  in  mind  that  the  surfaces  ranked  from  the  most  complex  to 
the  least  complex  (as  shown  by  the  power  spectra  given  in  Table  1): 

Surface  2 
Surface  1 
Surface  4 
Surface  3 
Surface  6 
Surface  5 

An  analysis  of  the  results  given  in  Tables  4  and  5  lead  us  to  the  following 
observations: 

1.  The  CISS  algorithm  performed  best  (according  to  the  standard  deviation  of 
the  error)  for  surface  5  and  worst  for  surface  4.  The  LISS  and  LIXY 
algorithms  performed  best  for  surface  5,  worst  for  surface  2  and,  for  the 
first  four  surfaces  listed  above,  the  algorithm  performed  as  the  CISS 
algorithm. 

The  ranking  of  the  surfaces  was  as  follows: 

CISS:  S,3,6,l,2,4 
LISS:  S,3,6,l,4,2 
LIXY:  5,3,6,1,4,2 

In  other  words,  all  three  algorithms  performed  best  for  the  three  least 
complex  surfaces. 


2. 


One  can  also  consider  each  surface  sei>arately  and  rank  the  algorithms 
according  to  how  they  performed  for  the  surface.  The  ranking  would  then 
be  as  follows: 


Surface  2: 
Surface  1: 
Surface  4: 
Surface  3: 
Surface  6: 
Surface  5: 


CISS, 

CISS, 

CISS, 

LIXY, 

CISS, 

LISS, 


LISS,  LIXY 
LIXY,  LISS 
LIXY,  LISS 
CISS,  LISS 
LIXY,  LISS 
LIXY,  CISS 


That  is,  for  surface  2  (the  most  complex),  CISS  performed  better  than  LISS 
which  performed  better  than  LIXY.  For  surface  5  (the  least  complex),  LISS 
performed  better  than  LIXY  which  performed  better  than  CISS.  To  put 
these  observations  in  perspective,  one  should  note  that  even  through  LIXY 
gave  the  worst  results  (compared  to  the  other  two  algorithms)  for  the  most 
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complex  surface,  it  may  have  produced  a  standard  deviation  of  the  error 
which  is  quite  acceptable  for  some  particular  purpose. 


3.  The  reason  why  the  three  algorithms  behaved  as  they  did  in  the  above 
observation  is  not  within  the  area  of  study  for  this  project.  (These  three 
algorithms  were  chosen  simply  because  they  were  available).  To 
investigate  the  reasons  why  the  algorithms  behaved'  as  they  did,  would 
involve  an  in-depth  investigation  of  the  relationships  between  the 
theoretical  basis  of  each  algorithm  and  the  physical  structure  of  each 
interpolation  region.  But  from  the  point  of  view  of  evaluating  algorithms 
in  a  practical  way  (by  using  synthetic  surfaces),  the  above  observation 
indicates  clearly  that  an  interpolating  algorithm  should  be  evaluated 
against  the  different  complexities  of  interpolation  regions. 


The  contour  interval  plays  a  significant  role  in  filtering  out  information 
contained  in  the  actual  terrain.  When  the  contour  interval  is  increased, 
information  about  the  original  surface  is  lost.  In  order  to  gain  an  insight  into 
this  asiiect  of  evaluation  of  interpolation  methods,  we  can  express  the  standard 
deviation  of  the  residuals  as  percentages  of  the  contour  interval.  This  is 
shown  in  Table  6. 


Table  6 

Standard  Deviation  as  a  Percentage  of  the  Contour  Interval 


Surface 

CISS 

LISS 

LIXY 

<T  (a) 

v/CIX 

»  (■) 

v/CIX 

9  (n) 

a/C  IX 

1 

70 

4.29 

20X 

9.64 

45X 

5.93 

28X 

2 

100 

7.22 

24X 

13.43 

44X 

14.15 

46X 

3 

20 

1.43 

23X 

1.66 

21% 

1.27 

21X 

4 

100 

7.99 

26X 

10.99 

36X 

10.28 

34X 

5 

10 

SOX 

0.67 

22X 

0.69 

23X 

6 

20 

1.44 

24X 

1.80 

SOX 

1.45 

24X 

An  examination  of  Table  6  shows  that  the  accuracies  of  the  three  inter¬ 
polation  methods  used  in  the  six  surfaces  range  from  20X  to  46X  of  the  contour 
interval.  Comparing  the  results  of  the  three  interpolation  methods  for  each 
surface  (that  is,  reading  along  the  rows  of  Table  6)  we  note  that  for  each 
surface,  one  of  the  methods  will  give  an  accuracy  of  between  20X  and  SOX  of 


the  contour  interval.  The  National  Map  Accuracy  Standards  give  a  limiting 
value  of  SOX  of  the  contour  interval.  Compared  to  these  standards,  all  three 
interpolation  methods  performed  well  on  all  six  surfaces. 

The  importance  of  this  research  is  shown  in  Table  6  and  in  the  comments 
made  in  this  section.  A  particular  algorithm  does  not  give  the  best  results  for 
all  types  of  interpolation  regions.  In  evaluating  algorithms,  an  important  factor 
(if  not  the  most  important  factor)  is  the  type  of  interpolation  region.  At 
present,  we  do  not  have  a  reliable  method  of  classifying  interpolation  regions. 
If  we  had  a  reliable  terrain  classification  method,  and  if  we  could  build  up 
experiences  of  using  different  interpolation  techniques  on  different  tyx>es  of 
terrain,  we  may  have  the  beginnings  of  an  Expert  System  in  this  field  of 
computer-assisted  mapping. 


6.2  Non-Numerical  Evaluation 

The  statistical  analysis  of  the  residuals  gives  a  measure  of  their  range  of 
magnitude  assuming  some  probability  level.  However,  what  matters  from  a 
cartographic  point  of  view  is  also  their  locations  and  their  pattern.  As 
mentioned  previously,  we  can  extract  information  concerning  the  existence  of 
blunders  or  any  kind  of  systematic  errors  and  draw  our  conclusions  for  the 
performance  of  the  interpolation.  For  this  reason  residual  maps  for  all 
datasets  and  for  all  the  interpolation  methods  have  been,  produced  and  are 
shown  below.  A  closer  look  at  these  maps  will  help  to  draw  more  general 
conclusions  about  the  performance  of  the  three  interpolation  methods. 

The  spatial  character  of  the  residual  maps  will  tell  us  something  about  the 
procedures  which  produced  the  residuals.  "Topography  type  terrain"  on  a 
residual  map  could  indicate  the  existence  of  systematic  error  in  the  process. 
Very  high  "peaks"  or  very  low  "depressions"  could  indicate  the  presence  of 
blunders.  Sparse,  low-valued  contours  could  indicate  the  random  nature  of  the 
process.  If  the  pattern  of  the  residuals  closely  resembles  the  pattern  of  the 
original  terrain,  one  would  expect  to  find  a  systematic  error  present. 

An  analysis  of  the  following  residual  maps,  leads  us  to  the  following 
general  observations: 

SURFACE  1:  The  best  results  are  given  by  the  CISS  residuals,  that  is,  we 
notice  a  general  random  pattern  of  low-valued  contours.  The  pattern  of  the 
residual  contours  does  not  resemble  the  contours  of  the  original  surface  -  this 
means  that  generally  there  is  no  systematic  error  in  the  interpolation.  We  do 
notice  three  areas  of  dense,  relatively  high-valued  contours.  These  areas 
appear  to  indicate  larger  errors.  Usually  blunders  would  be  flagged  out  by 
the  algorithm.  The  residuals  for  LISS  and  LIXY  show  larger  errors  indicating 
that  these  two  algorithms  give  worse  results  than  CISS.  This  is  an  expected 
result  -  linear  interpolation  is  not  suitable  for  complex  terrain. 

SURFACE  2:  The  residual  (error)  pattern  for  all  three  algorithms  appears  to 
be  random.  The  values  for  LISS  and  LIXY  are  leu*ge.  (Refer  also  to  Table  5 
where  the  maximum  absolute  errors  are  given).  The  beat  algorithm  is  again 
CISS  -  an  expected  result  as  Surface  2  is  a  complex  surface.  None  of  the 


altforithms  seem  to  have  any  ayatematic  errora  in  thia  aurface.  Again  we 
conclude  that  linear  interpolation  ia  not  auitable  for  thia  type  of  interpolation 
region. 

SURFACE  3:  Thia  aurface  ranked  fourth  in  order  of  complexity  of  the  aix 
aurfacea.  The  reaidual  mapa  indicate  that  LIXY  performa  the  beat.  None  of  the 
three  algorithma  ahow  evidence  of  ayatematic  error  -  that  ia,  the  contour 
patterna  are  of  a  random  nature.  Both  the  CISS  and  the  LISS  patterna 
indicate  a  large  error  in  about  the  aame  poaition  (the  denae  contoura  in  the 
lower  left  of  each  figure).  The  aize  of  the  error  ia  about  20  m  in  each  caae, 
approximately  equal  to  the  contour  interval  of  the  original  aurface.  Normally 
the  algorithma  would  flag  theae  areaa  for  further  inveatigation.  Our  concern 
in  thia  project  ia  not  with  the  actual  reaaon  for  theae  diacrepanciea,  we  are 
rather  drawing  attention  to  the  method  for  evaluating.  An  important  point  to 
note  ia  thia:  Table  5  indicatea  that  the  maximum  abaolute  error  for  CISS  for 
aurface  3  ia  20.43  m  and  the  reaidual  map  ahowa  that  thia  aize  of  error 
occurred  in  only  one  place  -  the  denae  contoura  in  the  lower  left  of  the  CISS 
reaidual  map  for  aurface  3.  Thia  fact  drawa  attention  to  the  value  of  having 
more  than  one  criterion  in  the  evaluation  proceaa. 

SURFACE  4:  Thia  aurface  waa  ranked  the  third  moat  complex  aurface.  The 
residual  mapa  show  that  all  three  algorithms  have  "pockets"  of  relatively  large 
errors,  but  CISS  gives  the  best  results.  Part  of  the  patterns  in  CISS  and 
LISS  are  the  same  and  generally  LIXY  has  a  different  pattern.  None  of  the 
residual  maps  resemble  the  original  synthetic  surface  -  an  indication  of  the 
absence  of  systematic  error. 

SURFACE  5:  This  surface  ia  ranked  as  the  least  complex  of  the  aix  synthetic 
surfaces.  One  would  expect  linear  interpolation  to  perform  well.  This  is 
confirmed  by  an  examination  of  the  LISS  and  LIXY  residual  maps.  Theae  maps 
illustrate  the  random  nature  of  the  errors  and  the  low  values  of  the  errora. 
On  the  basis  of  the  three  residual  maps  for  this  surface  one  can  see  that  both 
LIXY  and  LISS  give  better  results  than  CISS.  But  it  is  not  possible  to  tell 
which  of  the  two  linear  interpolatran  algorithms  is  the  better.  The  numerical 
values  given  in  Table  5  ahow  that  LISS  is  slightly  better  when  considering  the 
standard  deviation. 

SURFACE  6:  This  surface  was  ranked  as  the  next  least  complex  surface  after 
surface  5.  The  residual  maps  show  the  CISS  and  LIXY  algorithms  as  having 
better  fits  than  the  LISS  algorithm.  Here,  as  in  surface  5,  we  notice  that  the 
least  sophisticated  algorithm  (LIXY)  performs  well  in  terrain  of  low  complexity. 

The  residual  maps  have  not  indicated  the  existence  of  any  systematic  error 
in  the  interpolation  procedures  (if  the  residual  maps  resembled  the  original 
surfaces  or  exhibited  "topographic  shapes,"  then  one  would  suspect  the 
existence  of  systemmatic  error).  Attention  was  drawn  to  the  manner  in  which 
blunders  could  be  discovered.  The  actual  blunders  apparent  in  the  residual 
maps  were  corrected  so  that  the  remaining  random  errors  could  be  statistically 
analyzed.  In  an  evaluation  procedure,  one  may  need  to  investigate  the  cause 
of  the  blunders.  In  the  three  algorithms  being  used  for  this  project,  the 
blunders  were  due  to  a  bad  determination  of  the  contours  by  the  raster 
processing.  For  these  small  data  sets  they  were  easily  corrected. 


Figure  29 

Surface  4  -  CISS  ResidualB 
Scale  1:20,000  C.l.  =  20  m 


Figure  36 

Surface  6  -  LISS 
Scale  1:20,000  C.I.  =  2  m 


In  the  following  figures,  we  have  shown,  for  each  surface  and  for  each 
interpolation  method,  the  true  contours  (obtained  through  the  synthetic  surface 
equation)  and  the  contours  obtained  from  the  interpolated  grid  values.  Those 

figures  therefore  illustrate  the  practical  significance  of  the  error  statistics 

discussed  previously. 

As  mentioned  before,  the  contouring  algorithm  was  the  GSPP  package 

(Sunkel  1982).  The  accuracy  of  the  contouring  algorithm  will  not  affect  the 
conclusions  drawn  from  the  contour  maps  because  the  same  method  of 

producing  the  contours  is  used  for  both  the  original  surface  and  the 

interpolated  surface.  On  the  following  maps,  we  illustrate  the  practical 
differences  (on  a  scale  of  1:20,000)  between  true  contours  (obtained  from 
synthetic  surfaces)  and  derived  contours  (obtained  via  an  interi>olation 

algorithm). 

Our  discussion  of  errors  so  far  has  dealt  with  vertical  errors,  that  is, 
errors  in  height.  The  following  maps  translate  these  vertical  errors  to 

horizontal  errors  at  map  scale  for  each  algorithm  and  for  each  type  of 

interpolation  surface. 


Figure  38 

Surface  1  -  Origin*!  nnntour  and  via  CISS 
Scale  1:20,000  C.I.  =  40  ft. 


Surface  3  -  Ori 
Scale  1:20 
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Scale  1:20,000 
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Surface  5  -  Orii 
Scale  1:20 


Concluaions 


7.1  Method  for  Evaluation 

Our  investigation  as  to  methods  for  evaluation  have  shown  that: 

A.  Any  particular  algorithm  performs  differently  in  different  types  of 
interpolation  regions. 

B.  In  the  evaluation  of  an  algorithm,  one  needs  to  look  at  more  than  one 
criterion.  Each  criterion  gives  us  a  different  insight  into  the  performance 
of  an  algorithm.  That  is,  if  a  particular  algorithm  performs  well  for  a 
certain  interpolation  region,  it  may  not  perform  well  for  another  type  of 
interpolation  region.  Our  investigations  indicate  that  the  important  factor 
is  the  type  of  terrain  where  the  interpolation  takes  place  and  not 
necessarily  the  method  of  interpolation. 

At  this  stage,  our  proposals  for  a  method  of  evaluating  techniques 

(algorithms)  for  contour-to-grid  interpolation  would  be  as  follows; 

1.  Select  a  set  of  synthetic  surfaces  representing  different  types  of 
interpolation  regions. 

2.  From  the  equations  of  these  surfaces,  generate  "digitized"  contours  at 
suitable  intervals.  "Suitable  intervals"  would  be  contour  intervals  which 
are  likely  to  be  found  on  the  base  maps  to  be  digitized. 

3.  Apply  the  technique  (algorithm)  to  be  tested  to  the  "digitized"  contours 
and  obtain  predicted  grid  elevations.  The  grid  size  chosen  will  be  a 
function  of  the  base  'map  scale,  the  contour  interval,  and  the  user’s 
requirements. 

4.  Generate,  from  the  synthetic  surface  equation,  true  elevations  at  the  grid 
intersections  of  the  grid  used  in  3  above. 


5.  Using  the  results  of  3  and  4  above,  generate  the  following  statistics  and 

maps: 

a.  average  error; 

b.  standard  deviation  of  an  error; 

c.  maximum  absolute  error; 

d.  a  residual  map; 

e.  a  map  (preferably  in  two  colors)  showing  the  true  and  the  predicted 
contours.  The  true  contours  are  those  obtained  from  the  elevations  in 
4  above  using  a  contouring  package  and  the  predicted  contours  are 
those  obtained  from  3  above  using  the  same  contouring  package.  (The 
true  contours  should  be  the  same  as  the  contours  in  2  above.  Choose 
a  contouring  package  that  gives  this  condition.  Time  did  not  allow  us 
to  check  GSPP  -  of  more  importance  to  us  was  the  evaluation 
method). 


To  the  above  statistics,  one  may  have  to  add  others  when  using  large  data 
sets.  For  example: 

f.  computer  storage  space  required. 

g.  time  taken  to  perform  interpolation. 

The  above  statistics  and  maps  must  now  be  examined  and  a  decision  made 
on  the  performance  of  the  interpolation  technique  (algorithm).  One  would  be 
looking  for  evidence  of  systematic  error,  for  large  error,  for  blunders,  for 
horizontal  (mapped)  error,  etc.  as  described  in  Section  6  above.  The  final 
decision  will  depend  on  how  one  lists  the  following  requirements  in  order  of 
priority: 

*  accuracy 

*  speed  of  technique 

*  throughput  speed  (that  is,  including  preparation  time  and  time  to  achieve 
production  of  required  output) 

*  storage  space  requirements 

7.2  Future  Work 

A.  Our  investigations  described  in  the  previous  sections  have  shown  that  we 
do  not  have  an  essential  tool  for  the  evaluation  of  contour-to-grid 
algorithms.  This  tool  is  a  set  of  representative  synthetic  surfaces  -  that 
is,  a  suite  of  synthetic  surfaces  which  represent  the  types  of  interpolation 
regions  one  would  encounter  in  practice.  These  synthetic  surfaces  should 
be  capable  of  being  defined,  in  such  a  manner,  that  any  other  sampled 
surface  can  be  shown  to  be  similar  to  one  of  the  suites  of  surfaces 
mentioned  above.  Interpolation  algorithms  should  be  tested  on  the  suite  of 
representative  synthetic  surfaces  (as  described  in  7  above).  In  this  way, 
one  could  predict  the  behaviour  of  an  interpolation  algorithm  on  a 
peu'ticular  sampled  surface.  We  used  a  Fourier  series  for  setting  up  our 
(admittedly  not  entirely  satisfactory)  synthetic  surfaces.  We  indicated  that 
there  are  other  ways  of  generating  synthetic  surfaces.  Before  one  can  say 
that  the  Fourier  Series  method  is  not  suitable,  much  more  investigation  is 
needed.  With  our  time  limitations,  the  Fourier  series  was  used  because  it 
was  available  and  it  also  supplied  us  with  other  information,  like  the  power 
spectra.  A  high  priority  for  future  work  should  be  the  generation  o:' 
suitable  synthetic  surfaces. 

B.  In  order  to  evaluate  an  interpolation  technique,  one  also  needs  to  compare 
the  technique  with  other  "standardized"  techniques.  Some  of  these 
standardized  techniques  could  be  the  ones  used  in  this  study.  For 
example,  an  interpolation  technique  being  evaluated  may  give  accurate 
results  for  a  particular  terrain  type  but  a  standardized  technique,  like 
linear  interpolation,  may  give  results  within  the  accuracy  requirements. 

C.  This  study  has  not  considered  large  data  sets,  time  and  amount  of  storage 
space  required.  These  factors  should  be  considered  in  future  work. 
Considering  large  data  sets  will  bring  up  the  question  of  breaklines  and 
areas  of  homogeneity.  That  is,  even  though  large  data  sets  have  to  be 


px^essed,  the  interpolation  takes  place  only  between  breaklines  in 
homogeneous  areas. 

D.  This  study  has  only  considered  aspects  of  contour-to-^rid  interpolation. 
But  one  must  eventually  consider  also  grid-to-contour  and  grid-to-grid 
interpolation.  The  concepts  used  in  this  study  may  be  found  to  be 
applicable  in  these  cases. 

E.  An  investigation  which  is  related  to  the  above  considerations,  is  the 
question  of  the  density  of  the  data  collection.  That  is,  if  one  wants  to 
achieve  a  certain  accuracy  (of  interpolated  point  values  or  of  plotted 
contours)  using  a  particular  technique,  at  what  intervals  should  the  data 
be  collected?  Or,  if  the  data  is  very  dense,  does  one  have  to  use  all  the 
data?  The  results  of  this  study,  indicate  that  one  may  be  able  to 
determine  that,  for  a  particular  accuracy,  only  some  of  the  data  (for 
example,  every  fourth  contour  line)  needs  to  be  used. 
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APPENDIX  A 


BIBLIOGRAPHY  OP  SELECTED  REFERENCES 


Not  all  the  following  references  have  been  cited.  They  provide  a 
background  to  the  research  described  in  this  report. 


Bath,  Markus  (1974) 

"Spectral  analysis  in  Geophysics",  Elsevier  Sci.  Publ.  Co.,  Amsterdam  1974. 

It  is  an  extensive  presentation  of  the  mathematical  background  of  the 
Fourier  analysis,  together  with  many  applications  in  the  geosciences.  A 
discussion  of  topics  like: 

-  power  spectra/ their  computation 

-  reliability  and  presentation  of  spectra 

-  windows  and  filters  etc. 
is  also  presented. 

This  book  includes  an  extensive  bibliography  of  applications  of  spectral 
analysis. 

Boyle,  A.R.  (1980) 

"Scan  digitization  of  cartographic  data",  in  H.  Freeman  and  G.  Pieroni,  "Map 
Data  Processing",  Academic  Press,  1980. 

Deals  with  a  method  of  scan  digitization  of  existing  maps  in  the  light  of 
recent  advantages  in  software  for  line  thinning  and  vectorization. 

Clarke,  A.L.  (1982) 

"Interpolation  of  grid  digital  elevation  models  from  digitized  contours",  The 
Ohio  State  University,  Department  of  Geodetic  Science,  M.Sc.  thesis,  August 
1982. 

The  objectives  of  this  study  were  to  determine  the  best  interpolation 
algorithm  for  generating  high  fidelity  OEM’s,  and  to  determine  the  effect  on 
DEM  fidelity  of  the  inclusion  of  supplementary  non-contour  data.  A  new 
contour  -  specific  algorithm  (using  cubic  spline  interpolation  in  the  direction 
.  of  the  steepest  slope)  is  derived  from  an  optimum  fidelity  conceptual  model. 

Clarke,  A.L.,  Gruen,  A.  and  Loon,  J.C.  (1982) 

"The  application  of  contour  data  for  generating  high  fidelity  grid  digital 
elevation  models",  Auto-Carto  5  Proceedings,  August  1982,  pp.  213-222. 

The  special  characteristics  of  contours  as  a  DEM  data  source,  and  some 
published  contour-specific  interpolation  algorithms,  are  reviewed  in  this  paper. 
Test  results  of  these  algorithms  are  presented,  employing  both  synthetic  and 
real  surface  data. 
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Clericit  B.  and  Kubik,  E.  (1974) 

"The  theoretical  accuracy  of  point  interpolation  on  topographic  surfaces",  Intr. 
Soc.  of  Photogrammetry,  Commission  III,  Munchen  1974,  Proceedings  of  the 
symposium,  pp.  179-187. 


This  report  documents  an  investigation  into  the  theoretical  accuracy  of 
point  interpolation  and  volume  computation  for  topographic  surfaces.  For  this 
purpose  the  type  and  roughness  of  the  terrain  are  described  by  a  covariance 
function  of  exponential  type.  Based  on  this  model,  the  accuracy  of  interi)olated 
points  and  of  volumes  is  derived  by  the  application  of  the  law  of  propagation 
of  errors. 


Crochiere,  R.E.  and  Rabiner,  L.R.  (1981) 

"Interpolation  and  Decimation  of  Digital  Signals  -  A  tutorial  review", 
Proceedings  of  the  IEEE,  Vol.  69,  No.  3,  March  1981,  pp.  300-331.  ; 


A  tutorial  overview  of  multirate  digital  signal  processing  as  applied  to 
systems  for  decimation  and  interpolation.  Design  techniques  for  the 
linear-time-invariant  components  of  these  systems  (the  digital  filter)  are 
discussed. 


Crus,  Jaime  (1983) 

"Experiences  with  altimeter  data  gridding".  The  Ohio  State  University, 
Department  of  Geodetic  Science  and  Surveying,  Report  1347,  June  1983. 


The  effect  of  grid  spacing,  number  of  configuration  of  data  points  is  usod 
to  predict  a  grid  value,  and  prediction  method  used  is  illustrated.  The 
prediction  methods  tested  include  least  squares  collocation,  Bjerhammaur's 
inverse  distance  weighting  method,  Hardy’s  multiquadratic  method  using  .a 
hyperboloid  kernel,  a  general  harmonic  kernel  method  using  Krarup’s  kernel, 
and  a  piecewise  biquadratic  surface  fitting  method. 


Dahlquist,  G.  and  Bjorck,  A.  (1974) 

(translated  by  Ned  Anderson),  "Numerical  methods",  Prentice-Hall,  Inc.,  New 
Jersey  1974. 


Basic  numerical  methods  include  the  following:  basic  concepts  in 

approximation;  approximation  of  functions  by  Least  Squares;  use  of  polynomials; 
use  of  orthogonal  polynomials;  use  of  spline  functions. 


Davis,  D.M.,  Downing,  J.A.  and  Zoraster,  S.  (1982) 

"Algorithms  for  digital  terrain  data  modeling".  U.S.  Army  Engineer  Topo- 
grpahic  Laboratories  Report  number  ETC-0302,  July  1982. 


Author’s  abstract:  "Algorithms  and  test  results  are  described  for  three 
different  general  operations  used  in  computerized  modeling  and  mapping.  The 
algorithms  are  used  for  transforming  strings  of  contour  data  into  digital 
elevaticn  models  (DEM),  transforming  DEM’s  into  a  set  of  contour  strings  for 
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display,  and  smoothing  OEM's  using  several  filtering  and  convolution 
techniques.  The  use  of  DEM  smoothing  operations  for  the  simplification  of 
contour  lines  is  examined  and  two  methods  of  controlling  filter  effects  for  this 
applications  are  presented.  As  a  separate  research  effort  two  interpolation 
algorithms  which  can  be  used  for  resampling  OEM's  are  compared.  Visual  and 
statistical  measures  of  performance  are  presented.” 


Davis,  John  C.  (1973) 

"Statistics  and  data  analysis  in  Geology”,  J.  Wiley  k  Sons,  Inc.,  New  York  1973. 

In  addition  to  basic  statistics,  matrix  algebra  there  is  a  description  of 
interpolation.  Least  squares  regression  analysis,  autocorrelation  and  cross¬ 
correlation  analysis,  double  Fourier  series,  moving  averages  and  Kriging,  and 
statistical  testa  of  Geological  data.  An  extensive  list  of  FORTRAN  programs  is 
also  given. 


Davis,  Phillip  J.  (1975) 

"Interpolation  and  approximation”,  Dover  Publ.  Inc.,  New  York  1975. 

A  description  of  the  theories  of  interpolation  and  approximation,  Hilbert 
spaces,  orthogonal  functions,  orthogonal  polynomials,  measures  of  degree  of 
approximation  are  presented  in  extent. 

DMA  (1979) 

Defense  Mapping  Agency,  Product  Specificaticms  for  digitized  elevation  data  for 
Firefinder.  Hydrographic/Topographic  Center,  Washington  D.C.,  First  Edition, 
April  1979. 


Bbner,  H.  (1980) 

"HiFi  -  a  minicomputer  program  package  for  height  interpolation  by  finite 
elements",  14^^  Congress  of  ISP,  Hamburg  1980,  pp.  202-215. 

A  general  minicomputer  program  package  for  height  interpolation  is 
presented.  It  is  written  in  FORTRAN  and  interpolates  a  digital  height  model 
from  arbitrarily  disturbed  reference  points  and  breaklines,  using  the  Finite 
Element  method. 


Bbner,  H.  and  Reiss,  P.  (1981) 

"Experiences  with  height  interpolation  by  finite  elements",  ASP-ACSM  Fall 
Technical  Meeting,  1981. 

After  a  brief  description  of  the  HiFi  i>ackage  the  paper  gives  a  review  of 
the  experiences  gained  with  height  interpolation  by  finite  elements.  Reference 
is  made  to  the  accuracy  of  DEM  interpolation  and  to  the  operational  use  of 
HiFi. 


Bntfelis,  Theo  (1983) 

"Analysis  of  sea  surface  topography  using  Seasat  altimeter  data".  The  Ohio 
State  University,  Department  of  Geodetic  Science  and  Surveying,  Report  *343, 
March  1983. 

A  summary  of  Fourier  transform  and  low  pass  digital  filters  in  two 
dimensions  is  given.  A  geoid  corresponding  to  sea  surface  heights  (derived 
from  adjusted  Seasat  altimeter  data)  is  computed  up  to  a  minimum  wavelength 
of  2000  km,  using  the  GEMLl  gravity  model.  A  comparison  of  filtering  the 
spherical  harmonics  and  use  of  digital  filters  is  given. 


Bren,  Kamil  (1980) 

"Spectral  analysis  of  GEOS-3  altimeter  data  and  frequency  domain  collocation". 
The  Ohio  State  University,  Department  of  Geodetic  Science  and  Surveying, 
Report  *297,  February  1980. 

The  mathematical  background  in  spectral  analysis  is  applied  to  geodetic 
applications  is  summarized.  Frequency  domain  least-squares  collocation 
techniques  are  also  introduced  and  applied  to  estimating  gravity  anomalies  from 
GBOS-3  altimeter  data.  Power  of  each  frequency  for  discrete  and  non-periodic 
data  is  computed. 


Fernando,  K.V.  and  Nicholson,  H.  (1982) 

"Two-dimensional  curve  fitting  and  prediction  using  spectral  analysis",  IEEE 
Proceedings,  vol.  129,  Part  D.,  No.  S,  September  1982,  pp.  145-150. 

A  curve-fitting  or  prediction  problem  for  two-dimensional  or  cyclic 
processes  is  defined  and  solved  using  spectral  analysis.  The  assumed 
statistical  model  is  structurally  similar  to  the  Karhunen-Loere  expansion,  and 
the  technique  can  be  implemented  using  the  singular-value  decomposition. 


Franke,  Richard  (1979) 

"A  critical  comparison  of  some  methods  for  interpolation  of  scattered  data." 
Report  number  NPS-53-79-003,  Naval  Postgraduate  School,  Monterey,  California, 
1979. 

A  comparison  of  29  methods  for  the  solution  of  the  scattered  data 
interpolation  problem.  The  test  surfaces  used  were  very  smooth  -  when 
compared  with  topographic  terrain  encountered  in  practice.  The  given 
(reference)  points  were  randomly  scattered. 


Frederiksen,  P.  (1980) 

"Terrain  analysis  and  accuracy  prediction  by  means  of  the  Fourier 
transformation",  14^^  Congress  of  ISP,  Hamburg  1980,  Com.  IV,  pp.  284-293. 

In  this  study,  profiles  are  measured  in  different  geological  terrain  types, 
and  their  power  spectra  are  used  for  a  statistical  description  of  the  different 
landscapes.  By  means  of  the  spectrum,  the  standard  deviation  between  the 


terrain  surface  and  a  digital  terrain  model  is  estimated. 
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Prederikseni  P.|  Jacobi,  O.  and  Kubik,  K.  (1984) 

"Modelling  and  Classifying  Terrain”,  invited  paper,  Commission  III,  XVth 
Congress  of  the  International  Society  for  Photogrammetry  and  Remote  Sensing, 
Rio  de  Janeiro,  1984. 

Author’s  abstract:  "Different  methods  for  terrain  description  and 

classsification  are  presented.  Descriptive  calssifications  are  known  from 
geology.  These  descriptions  have  also  been  formalized  using  the  concept  of 
stochastic  processes.  Other  classification  methods  are  based  on 
Fourier-transformation,  classifying  the  terrain  types  according  to  their 
frequencies  of  undulation.  Recent  approaches  are  based  on  the  concept  of  self 
similarity,  observing  the  presence  of  lack  of  similarity  of  micro-  and 
macro-structures  in  the  terrain.  A  comparison  of  these  classification  models  is 
given  and  the  applicability  of  these  models  for  interpolation  and  other 
deductions  is  discussed." 


Hamilton,  Walter  C.  (1964) 

"Statistics  in  physical  science",  The  Ronald  Press  Co.,  New  York  1964. 

A  book  on  estimation  theory,  hypothesis  testing  and  least  squares  methods 
as  applied  to  data  of  physical  sciences. 

Hannah,  M.J.  (1981) 

"Error  detection  and  correction  in  Digital  Terrain  Models",  Photogram.  Eng.  & 
Remote  Sensing,  Vol.  47,  No.  1,  pp.  63-69,  1981. 

In  this  paper  algorithms  have  been  developed  to  detect  and  correct  errors 
in  digital  terrain  models.  These  algorithms  focus  on  the  use  of  constraints  on 
both  the  allowable  slope  and  the  allowable  change  in  slope  in  local  areas 
around  each  point.  Relaxation-like  techniques  are  employed  in  the  iteration  of 
the  detection  and  correction  phases  to  obtain  best  results. 

Hardy,  R.L.  (1977) 

"Least  squares  prediction",  Photogr.  Eng.  A  Remote  Sensing,  Vol.  3,  No.  4,  April 
1977,  pp.  475-492. 

The  basic  principles  of  least  squares  prediction  using  both  multiquadratic 
functions  and  covariance  functions  are  covered.  Multiquadratic  kernels  are 
based  on  geometric  and  physical  considerations,  rather  than  stochastic 
processes  as  in  the  case  of  covariance  kernels. 

Jacobi,  O.  (1980) 

"Digital  Terrain  model,  point  density,  accuracy  of  measurements,  type  of 
terrain,  and  surveying  expenses",  14^^  Congress  of  ISP,  Hamburg  1980,  Com. 
IV,  pp.  361-366. 


An  attempt  to  calculate  the  connection  between  the  accuracy  of  the  digital 
terrain  model,  the  terrain,  and  the  surveying  strategy,  by  means  of  a  sto¬ 
chastic  model  of  the  terrain  and  an  economic  model  of  the  surveying  expense. 


Jacobi,  Ole  and  Kubic,  Kurt  (1982) 

"New  concepts  for  digital  terrain  models".  Intern.  Society  of  Photogrammetry 
Commission  III,  Proceedings  of  the  symposium,  Helsinki  1982,  pp.  247-265. 

A  conceptual  model  for  the  toi>ographic  surface  is  presented.  A 
generalization  of  the  above  model  including  fractional  noise  is  attempted, 
together  with  a  small  discussion  on  the  concept  of  fractional  noise. 


Jancaitis,  J.R.  and  Junkins,  J.L.  (1983) 

"Modelling  irregular  surfaces",  Photogrammetric  Engineering,  Vol.  39,  No.  4, 
April  1973. 

A  two-independent  variable  interpolation  has  been  developed  for  modelling 
irregular  functions  of  two  variables.  This  method  represents  the  surface  as  a 
family  of  locally  valid  mathematical  functions  which  join  together  continuously. 


Kay,  S.M.  and  Marple,  S.L.  (1981) 

"Spectrum  Analysis  -  A  modern  perspective".  Proceedings  of  the  IEEE,  Vol.  69, 
No.  11,  November  1981,  pp.  1380-1419. 

A  summary  of  many  of  the  new  techniques  developed  in  the  last  two 
decades  for  spectrum  analysis  of  discrete  time  series  is  presented  in  this 
tutorial.  Techniques  discussed  include  the  classical  periodogram,  classical 
Blackman-Tukey,  autoregressive  (maximum  entropy),  moving  average,  maximum 
likelihood,  Prony,  and  Pisarenko  methods. 


Kearsley,  W.  (1977) 

"The  prediction  and  mapping  of  geoid  undulations  from  GEOS-3  altimetry".  The 
Ohio  State  University,  Department  of  Geodetic  Science,  Report  1267,  December 
1977. 

Geoid  heights  are  obtained  for  grid  points  in  the  region  to  be  mapped,  and 
two  of  the  parameters  critical  to  the  production  of  an  accurate  map  are 
investigated.  These  are: 

(i)  The  spacing  of  the  grid,  which  must  be  related  to  the  half-wavelength 
of  the  altimeter  signal  whose  amplitude  is  the  desired  accuracy  of  the 
contour,  and 

(ii)  the  method  adopted  to  predict  the  grid  values. 


Katsambalos,  Kostas  (1980) 

"Comparison  of  some  undulation  prediction  techniques  from  altimeter  data".  The 
Ohio  State  University,  Department  of  Geodetic  Science  and  Surveying,  Report 
•303,  July  1980. 
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Adjusted  GEOS-3  altimeter  data  are  used  to  predict  and  compare  geoid 
undulations  derived  from  two  collocation  models,  the  Bjerhammar  estimator,  and 
Hardy’s  multiquadratic  formula. 


Kraiky,  V.  (1977) 

"Grid-structured  reflexive  prediction  and  digital  terrain  modelling",  Amer.  Soc. 
of  Phot.,  Fall  Convention,  Little  Rock  1977,  pp.  447-454. 

If  the  prediction  process  is  based  on  a  covariance  function  separable  in 
coordinates,  conditions  are  set  for  an  effective  use  of  Rauhala’s  array  algebra. 
Its  advantages  are  demonstrated  for  the  "reflexive"  type  of  prediction 
implemented  with  or  without  least  squares  filtering. 

Kraiky,  V.  (1978) 

"DTM  interpolation  with  gliding  vectors".  Symposium  of  Digital  Terrain  Models, 
St.  Louis,  May  1978.  ^ 

Expanding  the  size  of  the  data  grid,  the  central  part  of  the  matrix  which  is 
inverted  in  the  process,  eventually  becomes  stabilized  at  certain  numerical 
values,  regardless  of  any  further  size  increase.  One  can  derive  a  table  of 
band  limited  vectors  to  serve  as  interpolation  or  smoothing  operators. 

Kraiky,  V.  (1980) 

"Spectral  analysis  of  interpolation",  14th  Congress  of  the  International  Society 
of  Photogrammetry,  Commission,  III  Hamburg  1980,  pp.  389-397. 

Basic  concepts  of  spectral  analysis  are  applied  to  interpret  interpolation, 
smoothing  and  parametric  transformation  based  on  uniform  sampling  as 
different  types  of  discrete  convolution.  Spectral  properties  of  interpolation 
are  then  discussed  with  the  main  emphasis  on  its  linear  least  squares  version. 


Lanczos,  Cornelius  (1966) 

"Discourse  on  Fourier  series",  Oliver  and  Boyd  Co.,  Edinburgh  1966. 

A  presentation  of  the  development  of  the  Fourier  series  as  well  as  the 
application  in  approximation  problems. 


Leberl,  F,  Kropatsch,  W.  and  Lipp,  V.  (1980) 

"Interpolation  of  raster  heights  from  digitized  contour  lines",  XIV  Congress  of 
the  International  Society  of  Photogrammetry,  Commission  III,  Hamburg  1980, 
International  Archives  of  Photogrammetry  Part  B3,  pp.  458-468. 

An  algorithm  for  converting  contour  information  to  digital  Elevation  Models 
is  given  here.  Examples  show  that  overall  interpolation  results  compare  well 
with  manual  interpolation. 
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Department  of  Geodetic  Science  and  Surveying,  The  Ohio  State  University, 
Columbus,  Ohio,  1982. 

An  introduction  to  OEMs  including: 

-  Data  acquisition 

-  Terrain  types 

-  Point  distribution  and  density 

-  Interpolation  procedures 

-  Overview  of  existing  Digital  Terrain  Models 

-  Applications  of  Digital  Surface  Models 


Moritz,  H.  (1978) 

"Introduction  to  interpolation  and  approximation",  in  H.  Moritz  and  H.  Sunkel 
eds.,  "Approximation  methods  in  Geodesy",  H.  Wichmann  Verlag,  Karlsruhe,  1978. 

The  following  topics  are  treated: 

-  Polynomial  interpolation  (formulas  of  Newton  and  Lagrange 

-  Least  Squares  approximation 

-  Polynomial  approximation  (Weierstrass  and  Chebyshev  theorems) 

-  Finite  elements 

-  Spline  functions 

-  Least  squares  interpolation  and  collocation 


Morrison,  J.L.  (1971) 

"Method-produced  error  in  Isarithmic  Mapping",  ACSM  Technical  Monograph,  No. 
CA-5. 

This  paper  deals  with  errors  introduced  by  isarithm  algorithms.  In  order 
to  evaluate  them  the  author  constructs  four  synthetic  surfaces  formulated  by 
logarithmic  functions,  trigonometric  series  or  equations  for  a  parabola.  Since 
these  surfaces  are  mathematically  defined,  an  absolute  error  measure  can  be 
defined. 


Rhind,  D.  (1975) 

"A  skeletal  overview  of  spatial  interpolation  techniques”.  Computer  Applications, 
Vol.  2,  No.  384,  pp.  293-309,  1975. 

A  critical  review  of  the  literature  concerning  interpolation  methods. 


Roberts,  P.A.  (1980) 

"Computer  based  technique  for  converting  a  contour  map  into  an  equispaced 
grid  of  points”,  U.S.  Dept,  of  Commerce,  National  Technical  Information  Service, 
September  1980. 

This  report  describes  a  technique  for  converting  a  contour  map  into  an 
equispaced  grid  of  points.  Included  are  discussions  of  the  problems  associated 
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with  the  digitization  of  a  contour  map  and  those  of  interpolation  in  general. 
Listings  of  all  programs,  which  are  written  in  computer  independent  ANSI 
FORTRAN  V,  are  given. 


Schultz,  H.H.  (1973) 

"Spline  analysis",  Prentice-Hall,  Inc.,  Englewood  Cliffs,  N.J.,  1973. 

This  paper  is  an  abbreviated  version  of  lectures  on  the  principle  and 
methods  of  approximating  functions  of  one  and  two  variables. 


Shut,  G.H.  (1974) 

"Evaluation  of  some  interpolation  methods".  Proceedings  of  the  symposium  of 
the  International  Society  for  Photogrammetry,  Commission  III,  Munchen  1974, 
pp.  160-170. 

Different  interpolation  methods  are  discussed.  For  example, 

-  Kraus’s  linear  least  squares  interpolation. 

-  Arthur's  interpolation  of  a  function  of  many  variables. 

-  Hardy's  interpolation  with  multiquadratic  equations 

-  Pointwise  interpolation 
and  so  on. 

Interpolation  methods  are  divided  into  six  groups. 


Shut,  G.H.  (1976) 

"Review  of  interpolation  methods  for  Digital  Terrain  Models",  the  CRnadian 
Surveyor,  Vol.  30,  No.  S,  December  1976,  pp.  389-411. 

A  review  of  the  methods  of  height  interpolation  for  digital  terrain  models 
that  have  been  published  in  photogrammetric  and  related  journals.  These 
methods  are  classified  into  six  groups,  according  to  the  type  of  surface 
representation  that  they  produce,  and  a  short  comparative  evaluation  is  made. 
In  addition,  the  construction  of  profiles  and  of  contour  lines  is  reviewed. 


Sunkel,  H.  (1980) 

"A  general  surface  representation  module  designed  for  Geodesy",  Report  No. 
292,  Department  of  Geodetic  Science  and  Surveying,  The  Ohio  State  University, 
Columbus,  Ohio,  1980. 

The  theoretical  background  and  the  manual  for  the  GSPP  package  designed 
for  interpolation,  least  squares  collocation,  generation  of  profiles  and  contours 
from  irregular  or  equally  spaced  reference  points. 


Tempfli,  K.,  and  Macarovic,  B.  (1978) 

"Transfer  function  of  interpolation  methods",  the  ITC  Journal  1978-1,  pp. 
50-79. 

The  objective  of  this  investigation  is  to  judge  the  performance  of 
interpolation  methods  in  general.  No  assumptions  are  made  with  regard  to  the 
type  of  terrain  relief,  but  sampling  is  assumed  to  be  homogeneous.  For  this 
purpose  a  special  analytical  approach  is  formulated  and  applied  for  the 
computation  of  Transfer  Functions. 


TempfU,  K.  (1980) 

"Spectral  analysis  of  terrain  relief  for  the  accuracy  estimation  of  digital  terrain 
models",  ITC  journal  special  issue,  (ITC  papers  presented  at  the  XIVth  ISP 
Congress,  Hamburg),  1980-3,  pp.  478-510. 

Deals  with  estimation  of  errors  introduced  in  profiles,  due  to  sampling 
interval,  interpolation  method,  aliasing  effect  etc.  The  theoretical  background 
of  spectral  analysis  used,  is  also  given.  Included  is  an  application  with  graphs 
in  logarithmic  scale. 


Werner,  H.  (1980) 

"Numerical  algorighms  for  interpolation  and  smoothing",  in  H.  Freeman  and  G. 
Pieroni,  "Map  Data  Processing",  Academic  Press,  1980. 


Wild,  E.  (1980) 

"Interpolation  with  weight-functions.  A  general  interpolation  method",  XIV 
Congress  of  the  International  Society  for  Photogrammetry,  Hamburg  1980, 
International  Archives  of  Photogrammetry,  part  B3,  pp.  780-793. 

A  discussion  on  the  use  of  weight  functions  in  interpolation.  It  includes: 

-  the  deduction  of  weight  functions  of  known  interpolation  methods. 

-  the  characteristics  of  weight  functions  and  their  effects  on  interpolation 
surfaces. 

-  weight  functions  as  a  linear  combination  of  base  functions. 

-  consideration  of  additional  information,  i.e.  curvature  lines,  lines  of 
steepest  slope,  highs,  lows,  etc. 


205  DIGITflLTERRRIME^Vftn^^ODEL  hUIlUUUUUUIU  Lflfi  w 
^  UMIV  RESEARCH  FOUNORTION  COLUMBUS  J  C  LOON  ET  RL 

JUL  85  0SURF-71E042  ETL-0292  DAAG29-81-O'0ie8 
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UNCLASSIFIED 


NL 


MICROCOPY  RESOLUTION  TEST  CHART 

1963-A 


ss>a«zaaaaassas3aa«ssssssssBma>ssaBSS3sssssssssssss=ss: 

I  C  I  S  S  ANDREW  L.  CLARKE  I 
I  CUBIC  INTERPOLATION  /  STEEPEST  SLOPE  /'  I 
I  /*  OllEAKLINE  DATA  HAY  BE  INCLUDED  /  SPRING  1902  I 


THIS  PROGRAM  USES  A  lUSTER  IMAGE  OF  CONTOUR  LINES 
'ro  INTERPOLATE  A  REGULAR  DEM.  THE  ALGORITHM  IS: 

"HEIGHT  AT  A  POINT  IS  FOUND  BY  CUBIC  POLYNOMIAL 
INTERPOLATION  WITHIN  FOUR  KNOWN  POINTS  ON  THE  I.INE 
OF  STEEPEST  SLOPE  PASSING  THROUGH  THE  UNKNOWN  POINT" 

STEEPEST  SLOPE: 

FOUR  DIRECTIONS.  ARE  TESTED  FOR  STEEPEST  SLOPE 
I.  WEST-EAST,  2.  NW-SE,  0.  NORTH-SOUTH,  4.  NE-SW. 

THE  SLOPE  IS  COHI'Ul'EO  FROM  THE  FIRST  CONTOUR  FOUND 
AT  EACH  END  OF  *I1IE  DIRECTION  LINE.  IF  THE  STEEPE.ST 
SLOPE  VALUE  OCCURS  IN  MORE  THAN  ONE  DIRECTION  LINE, 
THE  DIRECTION  HAVING  THE  SHORTEST  DISTANCE  BETWEEN 
CONTOURS  IS  USED. 

INTERPOLATION: 

AFTER  FINDING  THE  DIRECTION  OF  STEEPEST  SLOPE.  THE 
LINE  IS  EXTENDED  TO  THE  NEJCT  CONTOUR  AT  EACH  END. 

'IIIE  UNKNOWN  HEIGHT  IS  THEN  FOUND  DY  FITTING  A  CUBIC 
POLYNOMIAL  TO  'HIE  4  POINTS,  WITH  TOE  INTERPOLATION 
POINT  IN  THE  CENTRAL  SEGMENT:  Cl— C2— INT— C3— C4. 

IN  CASES  WHERE  EITHER  Cl  OR  C4  CANNOT  BE  FOUND,  AS 
AT  BREAK  LINES  OR  TOE  MAP  EDGE  ,  ONE  OF  TWO  INTER¬ 
POLATION  METHODS  MAY  BE  USED  (3  POINT  CASES): 

1.  FIND  Ce  OR  CS  AND  INTERPOLATE  BY  CUBIC  SPLINL. 

2.  FIND  ce  OR  CS,  THEN  CREATE  AN  EXTRA  POINT  AT 
EACH  END  BY  LINEAR  EXTRAPOLATION.  INTERPOLATE  DY 
CUBIC  SPLINE,  AS  IN  TYPE  1. 

EDGE  POINTS: 

INTERPOLATION  AT  TOE  EDGE  POINTS  IS  BY  CUBIC  POLY¬ 
NOMIAL  OR  SPLINE,  AS  ABOVE,  BUT  IS  LIMITED  TO  TOE 
DIRECTION  OF  TOE  EDGE. 

RASTER  DATA; 

THE  RASTER  IMAGE  IS  STORED  ON  DISK  (UNIT  «l),  WITH 
EACH  COLUMN  AN  UNFORMATTED  LOGICAL  RECORD.  EACH 
CORNER  HAS  A  HEIGHT  VALUE,  AND  THERE  ARE  NO  CAPS  IN 
THE  CONTOURS  IN  TOE  FOUR  TEST  DIRECTIONS.  IF  SPOT 
HEIGirr  DATA  ARE  INCLUDED.  THE  PIXELS  ALONG  TOE  LINE 
OF  STEEPEST  SLOPE  THROUGH  TOE  SPOT  SHOULD  CONTAIN 
INTERPOLATED  HEIGHTS. 

BREAK  LINE  DATA: 

IF  BREAK  LINES  ARE  INCLUDED,  THEY  SHOULD  HE  FLAGGED 
DY  HAVING  THE  CONSTANT  "DFLAC"  ADDED  TO  TOE  HEIGHT 
VALUE  AT  EACH  BREAK  LINE  PIXEL.  DATA  POINTS  ON  THE 
OPPOSITE  SIDE  OF  A  BREAKLINE  TO  A  DEM  POINT  WILL 
NOT  HE  USED  DURING  INTEIU’OLATION  AT  THAT  DEM  POINT. 

OUTPUT  DEM: 

THE  INTEIU’OLATED  DEN  IS  STORED  ON  DISK  UNIT  2,  Wmi 
EACH  COLUMN  AN  UNFORHAl'l'ED  LOGICAL  RECORD. 


'  •  T 
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INPUT  VARIABLES: 

THE  FOLLOWING  VARIABLES  MUST  BE  SET  BY  THE  USER  IN 
THE  NAIN  ritOGRAH.  DIMENSION  VARIABLES  MUST  AGREE 
WITH  THE  ACTUAL  DIMENSIONS  IN  THE  MAIN  PROGRAM. 


DIMENSIONS  OF  THE  INPUT  RASTER 
CONTOUR  INTERVAL  OF  THE  INPUT  DATA 
MAP  RESOLUTION  OF  THE  RASTER  DATA  (MM) 
CONTOUR  VALUE  CORRESPONDING  TO  MENU  NAV« 1 
GROUND  RESOLUTION  OF  THE  RASTER  (UNITS^CI) 


YES 


M.N 
Cl 
HES 
BASE 
RGD 
I BREAK 
BFLAG 
K.L 
I  DM 
TITLE 

suniioui' 

EDGE 
INSIDE 
DATA 
POLY 
SPLINE 
LINEAR 
TEST 
PRINT 
READ 
WRITE 

NOTES  i 

1.  THE  VALUE  ICRID»(  M- D/I  K- 1 )  SHOULD  BE  AN  INTEGER 
NUMBER.  THIS  IS  A  RESTRICTION  TO  THE  DIMENSIONS 
OF  THE  PRODUCED  DEN.  GIVEN  THE  DIMENSION  OF  HIE 
RASTER. 

2.  THE  CORNERS  OF  THE  INPUT  RASTER  SHOU'.O  HAVE  AT 
LEAST  APPROXIMATE  VALUES  (  IE. OTHER  THAN  ZERO) . 

TO  MODIFY  THEREFORE  THE  RASTER  VALUES  AT  HIE 
CORNERS  YOU  SHOULD  EITHER  INPUT  THESE  VALUES  IN 
THE  PROGRAM  BEFORE  THEY  WILL  BE  USED,  OR  DE’ITFIR 
YOU  CAN  USE  THE  PROGRAM  ■RASEDIT*  FOR  THIS  PUIU'OSE 


>  FLAG  FOR  BRFJIKLINE  DATA:  0  >  NO,  1  > 

>  CONSTANT  ADDED  TO  BREAK  LINE  HEICirPS 

>  dih»;nsions  of  the  output  den 

>  INTERPOIJVTION  MODEL  FOR  0  POINT  CASE 
:  60  CHARACTER  TITLE  READ  FROM  UNIT  5 

inks: 

-  CONTROL  ROUTINE  FOR  EDGE  INTERPOLATION 

-  CONTROL  ROUTINE  FOR  INSIDE  INTERPOLATION 
-FINDS  INl'ERPOLATION  CASE  AND  DATA  POINTS 

-  INTERPOLATES  WITH  CUBIC  POLYNOMIAL 

-  INTEIU'OLATES  WITH  NATURAL  CUBIC  SPLINE 

-  PERFORMS  LINEAR  INTERPOLATION 

-  COMPARES  INTERPOLATED  HEICRT  WITH  CONTOURS 

-  PRODUCE  A  LINE  PRINTER  HAP  OF  THE  RASTER 

-  READ  A  RASTER  OR  DEN  COLUMN  FROM  DISK 

-  WRITE  A  RASTER  OR  DEM  COLUMN  TO  DISK 


«  s  ss  s  8  s : 


SSSSSSSSSSSSSSS3SSSSSSSSrSSSS8SSSSSS8S8SSSSSSS 


I  L  I  S  S  AHDBEV  L.  CLARKE  I 

I  LIIIEAR  IlTTERPOLATIOir  /  STEEPEST  SLOPE  VIRTER  1982  I 

I  THIS  PROGRAM  OSES  A  RASTER  IMAGE  OF  CORTOOR  LINES  I 
I  TO  INTERPOLATE  A  REGIRJIR  DEM.  THE  ALGORITHM  USED  IS:  I 
I  HEIGHT  AT  A  POINT  IS  FOUND  BY  LINEAR  INTERPOLATION  I 
I  ALONG  THE  LINE  OF  STEEPEST  SLOPE  PASSING  THROUGH  I 
I  THE  POINT.  FOUR  DIRECTIONS  ARE  TESTED  FOR  STEEPEST  I 
I  SLOPE:  WEST  -  EAST.  NW  -  SE,  NORTH  -  SOOTH,  ME  -  SW.  I 
I  THE  END  POINTS  OF  THE  SLOPE  LINE  ARE  THE  FIRST  CON-  I 
I  TOUR  OR  EDGE  FOUND  IN  THE  TEST  DIRECTION.  I 

I  RASTER  DATA:  I 

I  THE  RASTER  IS  STORED  ON  DISK  UNIT  1.  I 

I  WITH  EACH  COLUMN  AN  UNFORMATTED  LOGICAL  RECORD.  I 
I  BACH  CORNER  HAS  A  HEIGHT  VALUE,  AMD  THERE  ARE  NO  I 
I  GAPS  IN  THE  CONTOURS  IN  THE  FOUR  TEST  DIRECTIONS.  I 

I  OUTPUT  DEM:  I 

I  THE  INTERPOLATED  DEM  IS  STORED  ON  DISC  UNIT  2,  WITH  I 
t  EACH  COLUMN  AN  UNFORMATTED  LOGICAL  RECORD.  I 

I  INPUT  VARIABLES:  I 

I  THE  FOLLOWING  VARl  ABLES  HOST  BE  SET  IN  THE  MAIN 
I  PROGRAM.  DIMENSIONS  MUST  BE  CONSISTENT.  t 

I  K.L  «  DIMENSIONS  W  THE  INTERPOLATED  DEM  I 

I  N.N  «  DIMENSIONS  OF  THE  INPUT  RASTER  I 

I  RES  >  RESOLUTION  (PIXEL  SIZE)  OF  RASTER  IN  HH.  I 

I  Cl  «  CONTOUR  INTERVAL  I 

I  BASE  >  CONTOUR  VALUE  ASSOCIATED  WITH  MENU  WAV  ■  1  I 

I  BFLAG  «  CONSTANT  ADDED  TO  BREAK  LINES  ( IF  USED)  I 

I  TITLE  a  «•  CHARAIHER  TITLE  READ  FROM  UNIT  5  I 

I  SUBROUTINES:  I 

I  READ  -  READ  A  RASTER  COLUMN  FROM  DISK  I 

I  WRITE  -  WRITE  A  DEM  COLUMN  TO  DISK  I 

I  EDGE  -  INTERPOLATE  EDGE  DEM  HEIGHTS  I 

I  INSIDE  -  INTERPOLATE  INTERIOR  DEM  HEI(»TS  I 

I  PRINT  -  PRODUCE  A  LIRE  PRINTER  HAP  OF  THE  RASTER  I 

I  NOTES  :  I 

I  l.THE  VALUE  ICRIDaCH-D/'lK-l)  SHOULD  BE  AN  INTEGER  I 
I  NUMBER.  THIS  IS  A  RESTRICTION  ON  THE  DIMENSIONS  I 

I  OF  THE  PRODUCED  DEM,  GIVEN  THE  DIMENSIONS  OF  THE  I 

I  RASTER.  I 

I  2. THE  CORNERS  OF  THE  INPUT  RASTER  SHOULD  HAVE  AT  I 
I  LEAST  APPROXIMATE  VALUES  (  IE. OTHER  THAN  ZERO)  .  I 

I  TO  MODIFY  THEREFORE  THE  RASTER  VALUES  AT  THE  I 

I  CORNERS  YOU  SHOULD  CITHER  INPUT  THESE  VALUES  IN  I 

I  THE  PROGRAM  BEFORE  THEY  WILL  BE  USED,  OR  BETTER  I 

I  YOU  CAN  I»C  THE  PROGRAM  *RA8EDIT'  FOR  THIS  I 

I  PURPOSE.  I 


L  I  X  Y  ANDREW  L.  CLARICE 

LINEAR  INTERPOLATION  /  X.Y  MEAN  WINTER  1982 

THIS  PROGRAM  OSES  A  RASTER  IMAGE  OF  CONTOUR  LINES 
TO  INTERPOLATE  A  REGULAR  DEM.  THE  ALGORITHM  USED  IS: 
HEIGHT  AT  A  POINT  HAY  BE  APPROXIMATED  BY  THE  HEAR 
or  THE  VALUES  FOUND  BY  LINEAR  INTERPOLATION  BETWEEN 
THE  NEAREST  CONTOURS  IN  THE  X  AND  Y  DIRECTION. 

RASTER  DATA: 

THE  RASTER  IS  STORED  ON  DISK  UNIT  1. 

WITS  EACH  COLUMN  AN  UNFORMATTED  LOGICAL  RECORD. 

EACH  CORNER  HAS  A  HEIGHT  VALUE.  AND  THERE  ARE  NO 
GAPS  IN  THE  X  OR  Y  DIRECTION  THROUGH  A  CONTOUR. 

OUTPUT  OEM: 

IHE  INTERPOLATED  DEM  IS  STORED  ON  DISK  UNIT  2,  WITH 
EACH  COLUMN  AN  UNFORMATTED  LOGICAL  RECORD. 

INPUT  VARIABLES: 

THE  FOLLOWING  VARI  ABLES  MUST  BE  SET  IN  THE  MAIN 
PROGRAM.  DIMERS lOra  HOST  BE  CONSISTENT. 

K,L  •  OIHERSIORS  OP  THE  INTERPOLATED  DEM 

M.R  >  DIMENSIONS  OP  THE  INPUT  RASTER 

RES  >  RESOLUTION  (PIXEL  SIZE)  OF  RASTER  IN  HH. 

Cl  «  CONTOUR  INTERVAL 

BASE  •  CONTOUR  VALUE  ASSOCIATED  WITH  MENU  NAV  >  1 
BFLAG  •  CONSTANT  ADDED  TO  BREAK  LINES  ( IF  USED) 
TITLE  •  «•  CHARACTER  TITLE  READ  FROM  UNIT  5 

SUBROUTINES: 

READ  -  READ  A  RASTER  COLUMN  FROM  DISK 

WRITE  -  WRITE  A  DEM  COLUMN  TO  DISK 

EDGE  -  CONTROL  ROUTINE  FOR  DEM  EDGE  INTERPOLATION 

INSIDE  -  CONTROL  ROUTINE  FOR  INTERIOR  INTERPOLATION 

LINT  >  LINEAR  INTERPOLATION  WITHIN  A  ROW  OR  COLUMN 

PRINT  >  PRODUCE  A  LINE  PRINTER  HAP  OP  THE  RASTER 

NOTES  : 

1. THE  VALUE  IGRIO«<M-l)/<K-l)  SHOULD  BE  AN  INTEGER 
NUMBER.  THIS  IS  A  RESTRICTION  ON  THE  DIMENSIONS 
OP  THE  PRODUCED  DEM,  GIVER  THE  DIMENSIONS  OF  THE 
RASTER. 

2.  THE  CORNERS  OF  THE  INPUT  RASTER  SHOULD  HAVE  AT 
LEAST  APPROXIMATE  VALUES  ( IE. OTHER  THAN  ZERO) . 

TO  MODIFY  THEREFORE  THE  RASTER  VALUES  AT  THE 
CORNERS  YOU  SHOULD  EITHER  INPUT  THESE  VALUES  IN 
THE  PROGRAM  BEFORE  THEY  WILL  BE  USED,  OR  BETTER 
YOU  CAN  USE  THE  PROGRAM  ’RASEDIT*  FOR  THIS 
PURPOSE. 


F  F  T  A  K  A  L 


PUllPOSF.  :  GIVEN  A  DATA  NATIIIX  OF  MXN  IT  COMPOTES  : 

1. HEAN  fi  VARIANCE  OF  THE  DATA 

2.  AVERAGE  POWER 

3.  POWER  PER  FREttUENCY  AND  ITS  %  CONTRIUUTION  TO  THE 
OVERALL  POWER 

4.  POWER  PER  HARMONIC  SURFACE 

5.  APPLIES  A  T-TEST  WITH  AX  SIGNIFICANCE  LEVEL  AND 
DELETES  THE  INSIGNIFICANT  FREOUENCIFiS 

6. CALULATES  (2).<3)  8  (4)  FOR  THE  FIL'IEREO  SURFACE 

7.  GIVES  BACK  THE  FILTERED  DATA 

DIMENSIONS  :  '/( H. N)  .Zl( M. N)  .ZSTARI M. N)  . XI M, N> 

P(Mt.NI)  .C(Ml.HI)  .PFIHl.Nl)  .CFini.Nl) 

IWK(  I)  .RWK(  I).  I«6»NAX(H.N)-i-15e 


NOTE  I  :  THE  USER  SHOULD  ALSO  CHANCE  THE  DIMENSIONS  OF  : 

X(H,N)  .CWK(M,N)  (COMPLEX  MATRICES)  IN  THE  ROUTINES 
"SPECTR*  a  ‘SPECTl" 

PARAMETERS  :  M  »  NUMBER  OF  ROWS  IN  DATA  MATRIX 
N  a  "  •  COLS  ■  *  • 

Nl«  NYQUIST  FREttUENCY  IN  X-AXIS 
N|a  •  ... 

THE  NVttUIST  FREttUENCIES  ARE  CALCULATED  BY  : 
Nl-N^2  IF  H*ENEN  NUMBER 

M|a(M+l)/2  :IF  MaODD  NUMBER 
A  a  SICIFICANCE  LEVEL  OF  T-TEST  (EG. 0.05) 

AUTHOR  :  I*.  C.  PAT  I  AS 
DATE  :  MAY  2tJ .  19H4 


SUnHOUTINF.S  USED  :  UTIL.SPECTR.SPECm  ,  FTEST,  HARM 
FFrjD.NDFD  (FROM  IMSL  LIBRARY) 


F  0  U  H  I  E  II 


PURPOSE  :  TO  CALCULATE  TUE  FOURIER  COEFFICIENTS  IN  THE  TRIGONO¬ 
METRIC  FOURIER  SERIES  EXPANSION  OF  A  SURFACE  GIVEN  RY 
A  HA'imX  OF  GRID  VALUES.  THE  COEFFICIENTS  WILL  COVER 
ONLY  UP  TO  3R0  HARMONIC  SURFACE 
ALGORITim  :  LEAST  SQUARES  ADJUSTMENT 

AUTHOR  s  P.C.PATIAS 

DATE  :  MAY6.I9U4 

DIMENSIONS:  THE  USER  NEED  ONLY  TO  CHANCE  THE  DIMENSIONS  OF  : 

A(NN.«7) .AT(49.NN) ,Z< NN) ,XX(NN) ,YY(NN) .AP(NN.49) 
WHERE  :  NN  «  TOTAL  NUMBER  OF  POINTS  OF  THE  SURFACE 
PARAMETERS:  THE  USER  SHOULD  CHANCE  THE  NEXT  PARAMETERS  : 

NN  «  TOTAL  NUMBER  OF  POINTS  OF  THE  SURFACE 
XO  *  X-COORD.  OF  THE  UPPER  LEFT  CORNER  OF  THE  GRID 
Y«  «  Y-COORD.  OF  THE  UPPER  LEFT  CORNER  OF  THE  GRID 
0  «  INCREMENT  DISTANCE  BETWEEN  THE  GRID  NODICS  ON 

X  a  Y  AXIS. 

OUTPUT  :  THE  EQUATION  FOR  THE  SURFACE  UP  TO  ORD  HARMONIC 
SURFACE  WILL  BE  GIVEN 

ROtrriNES  REQUIRED  :  THE  MATRIX  ROUTINES  OF  HIE 

SGEOlISC  I .  WYLEXEC'CSWMR  L I BRARY 


TERRAIN 


THIS  PROCRAH  CALCULATES  THE  Z-VALUES,  CIVEH  THE  FOURIER 
COEFFICIENTS  (UP  TO  3RD  HARMONIC  POWER)  AS  THEY  ARE  CONSTRU¬ 
CTED  BY  THE  PROGRAM  ■FOURIER*. 

IF  OPTION*  1  IS  CHOSEN  IT  ALSO  GIVES  A  PRINTER  CONTOUR  MAP 

DIMENSIONS  :  THE  USER  NEEDS  ONLY  TO  CHANGE  THE  DIMENSIONS  OF 
Z(M.N)  .W<  M.N) 

OPTIONS  :  0  *  GENERATION  OF  FULL  MATRIX 

IF  OPTION*e  IS  CHOSEN  SET  THE  UNNECl>!SARY 
PARAMETERS  TO  ZERO. 

I  =  GENERATION  OF  SPARSE  RASTER  CONTAINING  ONLY 
CONI'OUR  POINTS 

PARAMETERS  :  M  =  NUMBER  OF  ROWS  OF  THE  DESIRED  OUTPUT  MATRIX 
N  =  “  •  COLUMNS  «  »  ... 

XO*  X-COORI).  OF  'niE  LOWER  LEFT  CORNER 
YO=  Y-  ■  •  •  •  II  H 

lOPT  »  OPTION  USED  (  lE.e, 1) 

XX* DISTANCE  INTERVAL  ON  X-AXIS 
YY*  *  •  •  Y  » 

CI=  CONTOUR  INTERVAL 

BASE  a  Z- VALUE  OF  FIRST  CONTOUR 

TOL  a  TOLERANCE  VALUE  FOR  OBTAINING  THE  CONTOUR 
PIXELS 

FORMAT  !  THE  USER  SHOULD  ALSO  CHANGE  THE  FORMAT  ■  I  “  FOR 
THE  PROPER  READING  OF  THE  INPUT  COEFICIENIS. 

AND  THE  FORMAT  "ll*  FOR  THE  PRINTING  OF  THE 
RASTER  MATRIX 
AUTHOR  r  P.G. PATIAS 

DATE  :  MAY  22.1904 

SUBROUTINES  USED  :  PRINT  i 
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APPENDIX  C 

EQUATIONS  FOR  THE  SIX  SYNTHETIC  SURFACES 


For  all  the  following  equations: 

0  <  X  <  1000 
0  <  Y  <  1000 


Surface  1 


7. 


I 


»  10377.716 

♦  27. OK.  *  COS<Ky}  +  9.689  *  SlUfiXY)  * 

-140.006  WC03(2irf)  +  78.367  «  SIJt(2:C/)  + 

-90.210  CCC(GrCV)  +  0.106  *  SIHISKV)  + 

-19.091  COO’irCi)  +  -24.610  3'AUiai)  + 


15.0~3 

COOdGI)  s: 

COG  (IQ* )  + 

24.057 

GllldCC)  S- 

COGdCY)  + 

27.070 

••• 

•a* 

COGdCO  St 

003  (2E/)  + 

-40 . 439 

•M 

•a* 

GlNdCI)  ss 

oo3(2rao  + 

-34.543 

CCS  (Ed  S' 

C03(3ICy)  + 

-16.037 

«•« 

GlNdCi)  S' 

oo3(3:cy)  -i- 

-10.290 

<•» 

C03(2ia{) 

+  43.399  St 

-39.250 

St 

003  (2EC) 

St  C03dCY)  + 

-30,02.5 

•t* 

Gin(2EO 

S'.  COGdCi)  + 

14.397 

*** 

♦a* 

C03(2Ed 

St  C0S(2IC7)  + 

-37.590 

«•» 

*<* 

GCIKOiaO 

St  COGCJCT)  + 

1 1 .453 

a*« 

(:(;3(2:ao 

S:C0G(3IC/)  + 

JUKi'ICl) 

St  00.3(310')  + 

-~3.975 

C03(GI0d 

+  5.764  S: 

30.093 

t;OJ(3EO 

St  OOSdC/)  + 

71 ,077 

GIIKCJ'Oi) 

ss  COGdCY)  + 

-1 .151 

(:oG(3:3C) 

St  C03(2Ed 

-2.616 

'a> 

Gin  (3130 

ss  COG  (OCT)  + 

-21 .069 

a‘« 

C03(3iCO 

St  003(310/)  1- 

-26 . 306 

GIIKOEO 

*  COSdJICY)  + 

43,903  COS(KX)  *  SIIT(KY)  + 
20.213  s  SindiX)  s  SirU'CY)  + 
-37.742  *  COOliaC)  »  S1H(2KY)  + 

-12.417  s  sirfdco  #  sirrioKY)  + 

-13.399  COOdffi)  St  S!n(3KY)  + 

-22.C.34  SliTdK)  OIIdSI^Y)  + 

siniszx)  + 

10.537  ts  C0S(2:CI)  «  SindCY)  + 

-10.503  St  siruoic:)  ss  siikky)  + 

-30.072  St  COS(2IOO  St  SIIKOICY)  + 
-117.541  St  sinisrco  s:  Siri(2IC7)  + 
2.070  St  C03(2T<0  s:  GITHCIG:’)  ♦ 
77.697  ss  r;iri'r!TO{)  :.s  SUltCrCV)  + 

sinicro:)  + 

25.533  »  COS(OICt)  ss  SIIKKY)  + 
4.310  St  sin(3!a{)  St  stimc/)  + 
63.695  St  303(3100  s:  GlIKOICf)  + 
40.067  St  GIIKCIOO  St  G1W(2IC/)  + 
-10.294  ts  COGCOlOt)  St  GHKOIC/)  + 
-23.041  ss  SirdOICO  s:  SliUOIT/) 


Surface  2 


7.  a  1©S94.61S 

+  -‘JO.  139  #  COSIKY)  + 


34.713  »  S1B(KY>  + 


-5G.7L«3 

-C4.0Gn 

7.U19 
-3 .  '1'40 


cos(2rcY)  +  na.<.03  r.  siRimof)  + 

CC<3(C:C7)  +  -37.026  »  SITKGKY)  + 

cesdco  ■>  G2.GC-3  SINCICO  + 

COSdCI)  COSdC/)  12.603  a  COS(KX)  *  SIBIKY)  + 

SiriddO  a  COGdC/)  +  -IC.IOO  a  SINdQO  a  OIUdCY)  + 


.  iy;A 

C03(KI) 

•-  CCStS'C/)  + 

5.600  a 

COS  (ICO  * 

GTU'.CICY)  + 

-‘1’ .  'VJ 1 

3 III  (ICO 

cos  <  210')  •!■ 

-3.190  a 

SINdCE)  a 

GIIKSIO) 

-14. 6! 6 

••• 

COGdCn 

cos  (CEO)  + 

4. 200  a 

COSdCI)  a 

siKioro)  + 

■3.. 31 1 

••• 

GIIIdGO 

:  GOS(GIO)  + 

-3.764  a 

SuTdCI)  a 

Slj'dOIO) 

-60.310 

COS  (2100 

+  -203.447  a 

sndsrco 

+ 

-4.706 

cos(2rco 

•a  CeSdO)  + 

23.211  a 

COS  (2100 

a  sindO)  + 

-2.036 

5Ilf(2'aD 

a  CCS  (.'CO  + 

47.317  a 

SIIJ(2IOO 

a  STiKio)  + 

-07.214 

u* 

C03<2:Q{) 

a  cos  (210)  •*• 

-9.047 

a  cos(2:a{> 

a  SI  11(210) 

♦ 

or»/> 

4*< 

3IW(2EC[) 

a  CCS  (210)  + 

-22.190 

a  S!?r(2!0{,) 

a  3111(210) 

+ 

C03(2'O{) 

aCO.S(OIO)  •!- 

-7.20C  a 

COS  (2100 

a  SKTCOIO)  ■* 

-23 .  V'';0 

•  •• 

iii  11(2:30 

a  cc:j  (OiO/)  + 

-J.G50 

a  SIH(2IOO 

a  GIIKOlO) 

Jk 

140. 163 

«■« 

co3(o:go 

149.970  a 

SIIKCIGO 

+ 

16.029 

•  •• 

C03(0JCO 

a  COS  dO)  + 

-21 .027  a 

CCS(OIOO 

a  SKI  CO)  + 

-10.219 

GIIKOICO 

a  C0S(:O)  + 

-17.314  a 

SIIKClOO' 

a  sr.ido)  ■*- 

-34 .  .06 1 

C03(3kCC) 

a  C0J(21CC)  + 

-19.610 

a  COS(OIO{) 

a  r.iTKsroo 

29 . 975 

snKoico 

a  cos  (210)  r 

-20.912 

a  GIHOIOO 

a  SI II (2107) 

•f. 

40.1 12 

coGiorco 

a  C03(0r0)  + 

-’0.044 

a  C03(3KJO 

a  .SliKOIO) 

+ 

- 1 . 607 

*•* 

GHKOECO 

a  CCS(OIO) 

-40.640 

a  SINOtOO 

a  SKKGIO) 

f 


Surface  3 


o 
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Surface  4 


Z  *  S28.873 

+  1«.SI3  «  COS(Ky)  +  -26.S22  *  Slll(XY)  + 

84.037  «!C0S(2KY)  +  17.311  *  SIir(2KY)  + 

-77.004  «  COSISICV)  +  -«!-8.ll5  *  SINISICY)  + 

1.424  #  GOSCrCf)  +  -3.273  «  SIN(KX)  + 

-11.140  »  COSdCO  #  COSdO')  +  0.215  »  COSdOf)  *  SIlf(KY)  + 

-6.516  #  SINdCO  »  COSdCY)  +  -0.111  «  SIN(KX)  *  SIN(KY)  + 

-31.641  C  COSdCd  C0S(2ICY)  +  -11.723  *  COS(KX)  *  SIM(2KY)  + 

31.422  s  SINdCt)  !•;  C0S(2!CV)  +  16.233  *  SIN(KX)  «  SIN(2Ky)  + 

17.230  «  COSdGO  S!  C0S(3ICY)  +  -20.872  «  COS(KX)  «  SIN(3KY)  + 

-0.317  ss  SINdCO  COSISICY)  +  -34.470  «  SlNdOf)  SIN(3KY)  + 

13.369  »  C0S(2:(X)  +  20.835  «  SlN(2iaO  + 

24.209  «  C0S(2:{X)  *  COSdCY)  +  21.262  #  C0S(2KX)  *  SINIKY)  + 

3.992  #  SIN(2KX)  «  COSdOO  +  0.931  *  SIN(2KX)  «  SlNdt/)  + 

-32.397  #  COS(2IQO  «  C0S<2ICY)  +  -34.015  *  C0S(2ICX)  «  SIN(2KY)  + 

44.598  «  Sirf(2!00  »  C0S(2ICY)  04.224  *  SINIOIQO  «  SINIOKY)  + 

25.222  «  C03(2iOO  ^.-sCOSOKY)  +  -10.340  *  008(2100  *  SIN(3KY)  + 

03.133  *  SIW(2;CO  «  C0S(3ICY)  +  -11.538  *  Gin(2KX)  *  SINOKY)  + 

-46.120  *  C0S(3KX)  +  -17.337  w  815(3100  + 

-6.023  »  COS(3iaO  «  COSdCY)  +  -10.069  *  C0S(3ICX)  «  SINdCY)  + 

-10.734  SIN(3:00  J:  COSdCY)  +  14.239  *  SIN(3IGO  sx  SINdCY)  + 

7.297  *  COG(n:CO  C0S(2ICX)  +  16.001  V  C03(3iaO  *  SIN(2KX)  + 


-7.6.50  »  3IN(3iCX)  C03(2!CY)  +  -2.567  SIN(3IOO  *  SIN(2rCY)  + 
15.364  x:  C03(3I0O  :s  C0S(3KY)  +  -0.095  :X  COSOICX)  :S  SINOKY)  + 


a 


11.695  *  SIIKOICO  *  COS(aiCY)  +  11.934  #  SINOIQC)  «  SIN(3ICY) 
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Surface  5 


Z  >  388. 6«S 

+  0.796  *  COSdCY)  +  1.748  *  Sin<KY)  + 

1.726  «C0S<2KV)  +  -3.877  *  SIIf(2KY)  + 

3.622  #  COS(aKY)  +  3.876  #  SINOKY) 

0.427  *  GOSIKX)  +  0.444  St  3IH(KX)  + 

1.G42  #  COSdOO  S!  COSCITf)  +  0.237  *  COS(KX)  *  SUflKY)  + 

0.194  s:  SINdOf)  #  COSdCY)  +  -0.253  *  SIN(KX)  »  SIH(KY)  + 

0.112  *  COSdOO  *  C0S(2ICY)  +  0.284  *  COS(KX)  *  SIN(2ICY)  + 

0.525  *  SIWdOO  *  C0S(2ICY)  +  1,057  »  SIR(KX)  *  SIH(21CY)  + 

1.007  «  COSdOO  s;  C0S(3?CY)  +  -0.030  #  COS(KX)  *  SIN(3KY)  + 

0.199  St  SINdOO  COSOICY)  +  -0.005  *  SIN(KX)  «  SIN(3KY)  + 

-1.421  St  COS<2iaO  ♦  0.273  «  SIN(2KX) 

-0.5S2  St  COS<2IOO  *  COSIKY)  +  -0.268  *  C0S(2KX)  *  SIR(KY)  + 

0.656  *  SIN(2KX)  *  COS(KX)  +  0.066  *  Sirr(2KX)  *  SIIdKY)  + 

0.799  St  COS(2iaO  S'!  C0S(2KY)  +  0.032  *  C0S(2KX)  *  S1H(2KY)  + 

0.397  St  Sin(2iaO,S!  C0S(2KY)  +  -1.966  #  Sirf(2KX)  *  SIR<2KY)  + 

-1.497  St  C03(2IOO  stC0S(3ICY)  +  -1.329  *  COS(2IQO  *  SINOKY)  + 

1.903  s:  SIN(2IOO  st  COSOKY)  +  1.911  *  SIW(2IOO  *  SINOICY)  + 

1,216  s:  COSOIOO  +  1.130  »  SIIdOKX)  + 

1.460  St  COS(3lOO  s:  COSdCY)  +  0.376  #  C0S(3KX)  *  SIRdCY)  + 

0,223  St  SIW(3lOO  St  COSdCY)  +  0.417  *  SIR(3KX)  *  SIN(KY)  + 

1  .375  St  COSOfCX)  *  C0S(2KX)  +  -0.018  *  COSOKX)  *  SIR(2KX)  + 

0.136  St  SINOIOO  C0S(2ICY)  +  1.019  *  SINOKX)  *  SIN(2ICY)  + 

0,616  s:  COSiatOO  st  C0S(3KY)  +  0.735  *  COSIOIOO  #  SIIf(OKY)  + 

0.332  St  SINOKX)  *  C0S(3ICY)  +  -0,467  *  SIHIOKX)  *  SINOKY) 


Surface  6 


Z  »  766 

+  'i.502:*  COS(KT)  +  -1.442  *  SIIflKY)  + 

-6.592  «!C0G(2ICY)  +  -1.237  *  Sl»<2rC/)  + 

-1S.G?0  S'.  COSIQICY)  +  -9.214  S'.  3111(3107)  + 

0.7-55  S!  COS(IC{)  +  -3.015  s:  SIKdOl)  + 

-2.592  s:  COOdCO  S!  COS(K7)  +  -3.029  St  (K>9(KX)  *  SI1»(K7)  ♦ 

1.375  S!  3IN(IC{)  S:  COSdCV)  +  3.747  «  SiridOf)  Sin(KY)  + 

-3.0'?9  COSdOO  •:!  C0S(2:r/)  +  -1.416  «  COS(ICl)  *  SIIUGKY)  + 

-4.-570  s:  SltTdCD  s:  COS(2:OV)  +■  2.792  #  SllJdffi)  S:  SIN'.GICy)  + 

-0.415  SC'S  COSdCI)  *  COS(OICY)  +  7.507  *  COSdOO  SIIKCrCV)  + 

0.124  -S!  JindCO  #  COSCOICY)  +  -0.331  SlUdCO  *  GIIKCIC/)  + 

-12.614  S'.  GOS(2iaO  +  -0.365  #  S1NC2ICX)  + 

0.453  S'  COS(2IOO  «  COSdCV)  +  -2.973  #  C0S(2!aJ)  «  SiridCY)  + 

2.422  S'  3IN<2KX)  a  COSdOO  +  1 . 13C  a  SIIT(2:00  a  5in(T'")  + 

2.271  a  COS(2lOO  a  C0S(2ICY)  +  -4.490  a  C03(2K30  a  3111(2107)  ♦ 

0.0.3  a  3IN(2KX)  a  C0S(2K7)  +  -24.557  a  GIIT(.?fC:)  a  Sin(2K7)  + 

10.250  C03(2iaO  s:C0S(31C7)  +  2.659  a  C03(2IG{)  a  3111(3107)  ♦ 

2.SK'  a  3irr(2raO)  a  COS(3I07)  +  -5.447  a  GIN(2raf)  S:  .3IIT(.3I07)  + 

3.767  a  C03(3KX)  +  8.413  a  SIN(G!GC)  + 

0.492  a  003(3100  «  003(107)  +  1.267  a  003(3rCO  a  SITTd07)  + 

-1.115  a  3Ifl(3IGO  S'-  003(ICf)  +  0.463  st  3111(3100  a  SIVI(K7)  + 

-6.104  a  003(3100  a  003(2100  +  6.732  a  003(3100  S:  301(2100  + 

0.250  a  3IN(3!GO  a  003(2107)  +  5.958  a  SliKSIGO  a  Sirr(2KY)  + 

7.659  a  003(3100  a  C0S(3IC7)  +  -4.520  a  COCUOIGO  a  SKKOIC/)  + 

3.579  a  3(N(81GO  a  OOS(OKY)  +  -2.002  *  SII'KOKX)  S:  SIN(3K7) 


